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NOMENCLATURE 


a 

axis(M) 

Cab 

cf 

£,(«) 

fx 

poly'n 
polyn 
R n 
S{x) 
sat(x , a) 
X T 
X~ l 

X 

X 

X 

x (m) 

X a 

Q ab 

Si 

u J abc 


right-handed orthonormal axis system a 
axis operator; see equation (4-5) 

3x3 right-handed direction cosine matrix, rotation matrix from o to a 

binomial coefficient, namely k\/[i\(k — z)!] 

elementary (Euler) rotation about axis q through angle a 

partial derivative of / with respect to x 

polynomial of degree n with constant coefficients 

polynomial of degree n with possibly variable coefficients 

n-tuples of real numbers 

cross-product operator; see equation (3-3) 

saturation at a for |x| > a 

transpose of matrix X 

inverse of matrix X 

time derivative of x 

estimate of x 

vector x 

m th time derivative of x 

column matrix of coordinates of x with respect to a 

4 — ¥ 

Euler angles of axis a relative to axis b 
unit column matrix with 1 in location i 

c -coordinates of angular velocity of a relative to b 


SF m (x ) 

SF m { 0) 
SF m (l) 
( x,a ) 


Scalar Forms 

scalar form of x to order m, namely ( x, x , x ^\ . . . , x ^ m ^ ) 

zero scalar form, ( 0, . . . , 0 ) 
unit scalar form, ( 1, 0, . . . , 0 ) 
abbreviation for ( x, a, 0, . . . , 0 ) 
value in the k th location, namely x^ 


VF m (v a ) 
VF m ( 0) 
VF m {6i ) 
VF^Vah 


MF m (X) 

MF m { 0) 
MF m (I) 


Vector Forms 

vector form of v in *a to order m, namely ( v a , v a , . . . , Va^ ) 
zero vector form, namely ( 0, . . . , 0 ) 
unit vector form, namely ( <5j, 0, . . . , 0 ) 

scalar form of the i th component, namely (u a j, v a i , . . . , ) 

Matrix Forms 

matrix form of X to order m, namely ( X, X, . . . , X ^ m>) ) 

zero matrix form, namely ( 0, . . . , 0 ) 
unit matrix form, namely ( 7, 0, . . . , 0 ) 


v 


MF% 


RF™ 

AF ™(<*ab) 

PF 2 


matrix form of rotation matrix C ab , namely ( C ab , C ab , cffl ) 
Rotational Forms 

rotational form of a relative to b to order m, namely 

( C ab i u aba i ^ aba » • • • > w aba ^ ) 

identity rotation form, namely ( /, 0, . . . , 0 ) 

<-> *-* 

Euler angle form with sequence q of a relative to b to order m 
Euler parameter form of a relative to b to order m 


C f 

f[SF™{x)} 
f-l[ S F m (x)) 
T l f 

ji-1 

tay[SF m (x),r] 

MFyRF 

MF^RF 

MFyAF q 

MF^AF q 

MFyPF 

MF<PF 

XYyCY 

XY<CY 

XY^SP 

XY -<SP 


Functions of Dynamic Forms 

product operator 

control algorithm for pure feedback systems 

scalar form of f(x), namely ( f(x ), f(x), f( m ^(x ) ) 

scalar form of the inverse function / -1 (x) °f fi x ) 
algorithm computing the output of a dynamic system 
algorithm computing the control of a dynamic system 

first m + 1 terms of the Taylor series of x(t + r) 
transformation taking matrix forms to rotational forms 
transformation taking rotational forms to matrix forms 
transformation taking matrix forms to Euler angle forms in sequence q 
transformation taking Euler angle forms in sequence q to matrix forms 
transformation taking matrix forms to Euler parameter forms 
transformation taking Euler parameter forms to matrix forms 
transformation from Cartesian to cylindrical coordinates 
transformation from cylindrical to Cartesian coordinates 
transformation from Cartesian to spherical coordinates 
transformation from spherical to Cartesian coordinates 
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SUMMARY 


The formalism of dynamic forms is developed as a means for organizing and systematizing the 
design of control systems. The formalism allows the designer to easily compute derivatives to various 
orders of large composite functions that occur in flight-control design. Such functions involve many 
function-of-a-function calls that may be nested to many levels. The component functions may be 
multiaxis, nonlinear, and they may include rotation transformations. 

A dynamic form is defined as a variable together with its time derivatives up to some fixed but 
arbitrary order. The variable may be a scalar, a vector, a matrix, a direction cosine matrix, Euler 
angles, or Euler parameters. Algorithms for standard elementary functions and operations of scalar 
dynamic forms are developed first. Then vector and matrix operations and transformations between 
parameterization of rotations are developed in the next level in the hierarchy. Commonly occurring 
algorithms in control-system design, including inversion of pure feedback systems, are developed in the 
third level. 

A large-angle, three-axis attitude servo and other examples are included to illustrate the effectiveness 
of the developed formalism. All algorithms have been implemented in FORTRAN code. Practical 
experience shows that the proposed formalism may significantly improve the productivity of the design 
and coding process. 


1 INTRODUCTION 

This report presents a new procedure and a collection of algorithms for the solution of several 
problems associated with the design of automatic control systems. Our paradigm will be aircraft flight 
control, but the methods apply in other domains such as spacecraft attitude control, robotics, and process 
control. For flight-control design purposes, an aircraft may often be adequately modeled as a rigid body 
with force and moment generation that depends on the state of the motion of the rigid body, the controls, 
and wind. The design of the corresponding fully automatic flight-control system with large operating 
envelopes may be difficult for several reasons: 

For all but very restricted small-angle maneuvers, both . the fact that nonlinearities are associated 
with rigid rotation and the fact that the space of rotations is not flat become significant. Rigid body 
attitude is typically represented by direction cosine matrices or Euler angles in some sequence, or, in 
the case of spacecraft, Euler parameters (ref. 1 ). The link between the time derivatives of these attitude 
variables and the angular velocity and its derivatives is nonlinear and may even become singular. Thus, 
rotations introduce nonlinearities and singularity into the state equation. In order to avoid singularities 
it may be desirable to change from one representation to another at points along a flight maneuver. For 
example, the usual yaw-pitch-roll sequence becomes singular (gimbal lock) for 90 degrees of pitch, and 
in the vicinity of this condition it may be desirable to change to the yaw-pitch-yaw sequence. Each 
such change in representation entails a corresponding change in the state equations. The control system 
must be designed to operate in the various state space representations (coordinatizations), and the switch 
from one coordinate system to another must be made smoothly. Smooth patching of coordinates requires 
various-order derivatives of the right-hand side of the state equation (system function). 


In many cases the nonlinearity of the force and moment generators may not be ignored, especially 
for powered-lift aircraft that have strong nonlinear and rather complex interaction between power and 
aerodynamics. A typical algorithm for the computation of the total force and moment acting on the 
aircraft may contain more than 4,000 lines of FORTRAN code (ref. 2). The input to the algorithm is at 
least 19-dimensional, consisting of the state, which is at least 12-dimensional, plus the controls, which 
are at least 4-dimensional, plus wind; the output is 6-dimensional, consisting of the 3-dimensional force 
and moment vectors. Inside the algorithm, the input flows through many successive functions so that the 
analytic form of the multivariable function represented by the algorithm is a deep nesting of elementaiy 
functions and table interpolations. The depth of nesting (e.g., square of sine of square root of sum of 
squares of . . . etc.) easily exceeds dozens of levels. Consequently the analytic computation of such a 
simple object as the overall 6 x 19 input-output Jacobian matrix may be a formidable task. But such 
mathematical objects are needed if the designer is to improve system performance by taking advantage 
of the information contained in the force and moment model. 

The number of controls often exceeds the basic four. In addition to the three moment and one 
throttle controls, there may be controls for directing thrust (such as a one- or two-degree-of-freedom 
nozzle), direct lift (such as a spoiler), side force, and flaps. The set of all the controls may be redundant 
in the sense that many combinations of controls produce the same total force and moment on the aircraft. 
This redundancy may be resolved advantageously by partitioning the set into two sets: the nonredundant 
set of active controls and the set of parametric controls. The active controls are manipulated by the 
regulator; the parametric controls are manipulated by a configuration-management system so as to 
maintain selected control margins for the active controls as well as to maintain certain functions of the 
state (such as the angle of attack) within the assigned limits (ref. 3). Each such partition represents 
a particular control mode. Typically there are several control modes, and there may also be several 
tracking modes. Each such mode is associated with a particular functional dependence of the output on 
the state. For example, near hover, the output may be the three Cartesian coordinates of the velocity 
vector; at a higher speed, the output may be defined as cylindrical coordinates of the velocity, namely, 
horizontal speed, heading, and vertical speed; at a still higher speed, the output selected for tracking 
may be the spherical coordinates of the velocity, namely, airspeed, glide-path angle, and heading angle. 
Each combination of control and tracking modes defines an operating mode. Thus, near hover, the 
nozzle may be an active control, and the pitch angle may be programmed independently. At high 
speed, the nozzle may be fixed, that is, programmed independently, and regulation is then achieved 
indirectly through the pitch angle. Each operating mode is a separate control problem with its particular 
control variables, output variables, and, possibly, state variables and state equation. The control -system 
design must incorporate many such operating modes and provide smooth intermode transitions. Smooth 
patching of modes requires various-order derivatives of the system function. 

The flight-control system as considered in this report includes the functions of configuration man- 
agement, in which operating modes and reconfiguration commands are computed; guidance, in which 
flyable reference trajectories linking way points given by, say, air traffic control, are generated; and 
regulation, which ensures tracking of reference trajectories in spite of unavoidable uncertainties and 
approximations in system modeling. If the operating envelope is small enough relative to system non- 
linearity, then linear design methods based on a single Jacobian matrix (perturbation model) of the 
system function, evaluated at a single operating point (trim point), are adequate for the purposes of 
regulator design. In such a case, the Jacobian matrix may be computed numerically by perturbing each 
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input variable (ref. 4). For larger envelopes, robust single-point linear designs (ref. 5) based on one 
representative value of the Jacobian matrix may be adequate in spite of actual variations of the matrix. 
Since the operating point is now a variable, feed-forward signals (particular solutions) may have to be 
provided to reduce tracking error and unload the feedback (ref. 6). Nonlinear methods become essential 
for the design of flight-control systems with large operating envelopes. Two situations arise: If the 
inner-loop dynamics (attitude control) can be made sufficiently fast relative to the outer-loop dynamics 
(trajectory control), then the nonlinearity of the force and moment function may be removed by means 
of numeric inversion (refs. 7-13). The nonlinear control theory based on differential geometry (refs. 14 
and 15) provides design techniques when such a separation is impossible. One fruitful technique is 
based on a coordinate change of state and control in order to simplify the form of the state equation. 
In certain practical cases, a coordinate change may suppress the nonlinearity to the extent that linear 
design techniques become applicable in the new coordinates (ref. 16). Techniques are also available 
for the generation of the nonlinear analogs of the particular solution (refs. 17 and 18). The practical 
drawback of such techniques, for the case of flight control, is that they have a voracious appetite for 
various-order derivatives of the force and moment functions. 

Thus, the design of large-envelope flight-control systems is difficult because the state space is 
not flat, the force and moment function is big and complicated, and many operating modes must be 
considered. The difficulty can be further traced to the need for high-order differentiation of the system 
function. There are three differentiation techniques: hand, symbolic, and automatic. Hand differentia- 
tion and coding is very tedious and highly unreliable for the size of problems being considered. The 
other two alternatives are much more appealing. The symbolic approach would be to machine translate 
the, say, FORTRAN code for the system function into the appropriate language (such as MACSYMA, 
MATHEMATICA, or MAPLE) within which differentiation is defined, and then proceed with the non- 
linear design techniques employing the derivatives. Applications of this approach to relatively small 
systems have been successful (ref. 19). However, for larger systems (4,000 lines of FORTRAN) involv- 
ing deeply nested functions, symbolic methods may be slow and may often produce large, unmanageable 
expressions (ref. 20). 

The remaining choice for the computation of derivatives, automatic differentiation, is based on the 
fact (known since Leibnitz) that Taylor series, which carry derivatives as coefficients, can be propagated 
through an arbitrary sequence of elementary functions without any truncation error. Thus, automatic 
differentiation does not suffer from the rapid-chain-rule fanout of terms that plagues the symbolic dif- 
ferentiation. Furthermore, machine translation into an automatic -differentiation language is as practical 
as it is for symbolic languages (ref. 21). 

The theory of dynamic forms (ref. 22) described in the present report may be considered to be a 
particular example of automatic differentiation and a basis for a formal language for the computer-aided 
design of automatic control systems. 

A dynamic form is defined as a variable together with its time derivatives up to some fixed, but 
arbitrary, order. The variable may be a scalar, a vector, a matrix, a direction cosine matrix, Euler angles, 
or Euler parameters. Most of the report is devoted to the translation of a set of elementary functions and 
operations into corresponding functions and operations on dynamic forms. The set is rich enough so 
that typical system functions occurring in flight control may be assembled from the members of this set. 
Whereas many examples are provided to demonstrate the application of the methodology to automatic 
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control, the emphasis in the present report is on the formalism of dynamic forms. The emphasis in 
future reports will be on the application of dynamic forms and on the reformulation of control problems 
in terms of dynamic forms. 


2 SCALAR FORMS 

Mathematical models of practical dynamic processes such as an aircraft frequently contain functions 
of functions to many levels. A small example of a typical case is shown in figure 2.1. In the figure the 
output y is related to the input x = (xi,X2,X3) by a function /, which is built up from the elementary 
functions such as trig functions, powers, roots, exponentials, and, of course, addition, subtraction, 
multiplication, and division. The function / so constructed then becomes a block in a higher level 
function and so on until the final level, say w = F(u), is reached. Implementation of such a function F 
as an algorithm on the computer is routine, even though in practice the algorithm may easily take 
4,000 lines of FORTRAN. On the other hand, the symbolic expression of F becomes unwieldy. Even 
the simple example fragment f of the complete algorithm F is beginning to look complicated. 

y = f(xi,x 2 ,x 3 ) = fc{h(x 3 ) * UlfiM F f 2 (x 2 )]}/fa(x 3 ) (2-1) 

The function is six levels deep in the sense that six function calls (/, fa, *, fa, +, f 2 ) are needed to get 
from x 2 to y. Now, whereas the analytic form is not needed for simulation, it is often used during 
design and analysis for the computation of, for example, time derivatives, gradients, partial derivatives, 
and Jacobian matrices. Suppose that we wish to compute the first five time derivatives of the output y 
in figure 2.1 given the input (x\,x 2 ,x 3 ) and its first five time derivatives. Repeated use of the chain 
rule would produce rather long expressions. In general the length grows rapidly with the depth of the 
nesting. The same considerations apply to partial and other derivatives since they can be expressed in 
terms of time derivatives. Indeed, an effective procedure for computing time derivatives can be easily 
adapted for computing other derivatives. This procedure will be discussed later in the report, so we 
focus on the computation of time derivatives. The scalar dynamic forms discussed next greatly simplify 
the computation of time derivatives. 

Suppose that a scalar variable x is a function of time. Imagine an array with x in location 0, time 
derivatives running to the right, and time integrals to the left: 

(..., ffx, Jx, x, x, x, x@\ ...) (2-2) 




x 2 


x 3 



Figure 2.1. Typical nested fragment /. 
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The scalar form SF m (x ) of order m is a glimpse at this array through a window at locations 0 through 
m. That is, 

SF m (x) = (x,x,...,xW) (2-3) 

where x € Z? 1 is a function of time and x( m ^ in location m denotes the ra^ time derivative of x. 
Location 0 contains x = x^. Scalar forms will be denoted by the symbol SF to distinguish them from 
other dynamic forms to be introduced later in the report. 

A time derivative of a scalar form will be defined as a right shift of the window: 

d/dtSF m (x ) = SF m (x ) = a;( m+1 ) ) (2-4) 

The integral of a scalar form is given by a left shift of the window: 

J SF m (x) = SF m (fx) = (fx,x,x,...,x( m - 1 )) (2-5) 

We use the following notation for extracting a time derivative from a form: 

xW if 0 < k < m 

{SF m (x)} k = j (2-6) 

„ 0 else 

There is a close relation between a dynamic form SF m (x ) and the first to + 1 terms of a Taylor series 
expansion of x(t) at t. Thus, 

x(t + r) « x + xt + + . . . -I -x^ m V m 

v 2 ml 

where the time derivatives are contained in SF m (x). We use the notation 

171 u 

x(t + r) * t] = £ {SF m (x)} k T k /kl (2-7) 

A;=0 

Algorithms such as equation (2-1) may be considered to be of order zero. The primary objective of 
this report is to develop a formalism for easy conversion of such zero-order algorithms into corresponding 
algorithms of order m > 0. Thus the algorithm in (2-1) will translate into 

SF“( y) = h (f 3 {SF m (:r 3 )] . /„{/i [Sf™^)] + / 2 [S/’ m (* 2 )]}) /APf-fo)] < 2 ' 8) 
where the input is given by the three scalar forms (variables and their time derivatives up to order m) 

SF m (x 2 ) = (x 2 ,± 2 ,...,x^ m) ) 

SF m (x 3) = (x 3 ,i3.-.-.4 m> ) 
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the output is given by the scalar form (variable and its time derivatives up to order m) 


SF m ( V ) = (y,y,. 

and the meaning of fj and the operators in equation (2-8) will be developed next. The formalism 
will allow us automatically to translate large zero-order algorithms such as an aircraft total-force-and- 
moment subroutine, which, as noted earlier, may contain 4,000 lines of FORTRAN, into an order 
algorithm. That in turn may be used for the design of control algorithms such as the generation of 
reference trajectories and model inversion. Examples of such designs will be given later in the report. 

Let us now proceed with the development of the formalism. We shall first convert elementary 
functions to the corresponding functions of scalar forms. Then, with these functions as basic building 
blocks, we will assemble a hierarchy of composite functions that are particularly useful for the design 
of control systems. 


Algebra 

Zero and unit scalar forms are defined, respectively, by 

(2-9) 
( 2 - 10 ) 


£F m (0) = (0, ...,0) 
SF m (l) = ( 1, 0, . . . ,0) 


The scalar form of time is 

SF m (t) = (t, 1, 0,..., 0) 
To simplify notation, we assume padding with zeros: 


(x,a) = (x, a, 0, ...,0) 


Thus, we may write 

SF m ( 0) = (0) 
SF m ( 1) = (1) 
SF m (t) = (t, 1) 


( 2 - 11 ) 

( 2 - 12 ) 


It is possible to define sum, product, inverse, and division for scalar forms. If z = ax + by and a 
and b are constant in time, then for 0 < k < m 

z W = ax^ + by ^ (2-13) 

An outline of the algorithm for scalar-form sum is shown next, where scalar forms are treated as scalar 
arrays. 
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ALGORITHM: SF m (z ) = SF m (x ) + SF m (y ) 
do k = 0, m 
z^) = a * x^) b * 
end do 


The effect of the algorithm will be denoted as 

SF m (z ) = aSF m (x ) + bSF m (y ) 


Note that 

SF m (x) - SF m (x) = SF m ( 0) 


(2-14) 


If z = xy, then the k th 
ref. 23): 


derivative, 0 < k < m, of z is given by the Leibnitz product rule (convolution. 


= £ Cfx^-VyW = £ CfxMyib-i) 

i—0 i=0 


(2-15) 


where the binomial coefficient Cf = k\/[i\{k — 


i)!] may be computed by means of the Pascal triangle: 
for k > 0 


(2-16) 


[ C\ = C\Z\ + Cf _1 for 2 < fc and 1 < i < k - 1 

An outline of the product algorithm is shown next, where as before scalar forms are treated as scalar 
arrays. 

ALGORITHM: SF m (z ) = SF m (x) * SF m (y ) 
do k = 0,m 
= 0 

do i = 0, k 

z^) — zW + Cj * * y(*) 

end do 
end do 


The effect of the algorithm will be denoted as 

SF m (z ) = SF m (x) * SF m (y) (2-17) 

The scalar form product commutes since it commutes for real numbers: 

SF m {x ) * SF m {y ) = SF m (y) * SF m (x ) 


It may be noted that the Leibnitz rule holds also for objects other than scalars. It matters only that 
the product is defined for objects x and y so that its time derivative 

(xy)( 1 ) = x^y -1- xt/ 1 ) (2-18) 

We will take advantage of this fact later when we consider vectors and matrices. For now, we return to 
scalars. 
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Powers and Polynomials 

The conversion of a zero-order algorithm to the corresponding algorithm of order m is illustrated 
by the following very simple case. 

If z — x n for an integer n > 1, then an obvious algorithm for computing z is given by 

ALGORITHM: Z = X n 
Z = 1 

do i = 1 , n 
z — z * x 
end do 

The rule for conversion to order m is simple: replace any variable of time by its scalar form. Thus, the 
algorithm for computing 

SF m (z) = [5F m (x)] n (2-19) 

is given by 

algorithm: SF m (z ) = [5'F m (a:)} n 
SF m (z) = SF m ( 1) 
do i = 1,71 

SF m (z) = SF m (z)*SF m (x) 

end do 

where the algorithm (subroutine) for SF m (z ) * SF m (x) has been already constructed (see eq. (2-17)). 
The algorithm works for any x but requires n scalar-form products. For x ^ 0 an algorithm requiring 
essentially one product will be constructed when fractional powers are considered later in the report. 

Polynomials occur frequently enough in practice to deserve consideration. Let poly' n (a,x ) denote 
a polynomial of degree n in x with constant coefficients a = (ao, . . . , a n ), that is, 

poly' n (a, x) = ao + a\x + . . . + a n x n = ao + (ai + . . . + (a n _i + a n x) ...x) (2-20) 

If 2 = poly' n (a , x), then the following algorithm is a possible realization of poly ' n : 

ALGORITHM: 2 = poly' n (a,X ) 

2 = a n x 
do i = 1, n — 1 
2 = (a n _j + z) * X 
end do 
z — ao + z 

The order of the algorithm is raised from zero to m simply by replacing 2 and x by SF m (z) and 
SF m (x), respectively: 
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algorithm: SF m (z ) — poly' n [a, S F m (x)] 

SF m (z) = a n SF m (x) 
do i — 1, n — 1 

SF m (z ) = [a n _jSF m (l) + SF m {z)} * SF m (x) 

end do 

SF m (z) = a 0 SF m (l ) + SF m (z ) 

This algorithm will be denoted as 


SF m (z ) = poly' n [a,SF m (x)\ 


( 2 - 21 ) 


ft 


Example 2.1. Suppose that we need to generate the following polynomial function of time: 

x = poly 4 (a, t) = l + t- t 2 1 2 + t 3 /3! + i 4 / 4! 
and that we need derivatives to order six. Then we need to compute 

5F 6 (x)=pol^[a,5F 6 (f)] 

Thus, for example, at t = 1, the scalar form of time is 

SF Q {t) = (1.00, 1.00, 0.00, 0.00, 0.00, 0.00, 0.00) 
and a call to the poly algorithm with degree= 4 and coefficients 

a= ( 1,1 ’”2’3!’4 t) 

produces the scalar form of x, 

SF 6 (x) =poly' 4 [a,SF 6 (t)] = (1.71, 0.67, 0.50, 2.00, 1.00, 0.00, 0.00) 

This result may be checked by hand: x(l) = 1.708 . . ., x(l) = (1 — t + 1 2 / 2 + f 3 /3!)t=i = 0.666 . . etc. Next 
suppose that we must pass x through a nonlinear block represented by a polynomial of degree three in x: 

z = polys(b, x) = 1 — 0.50x — 0.50x 2 + O.lx 3 

Then another call to the poly algorithm of degree three and coefficients 

b= (1, -0.5, -0.5, 0.1) 


gives the scalar form of z, 

SF 6 (z)=poly' 3 [b,SF 6 (x)] = (- 0.82, -0.89, -0.66, -2.46, -0.38, 7.17, 30.31) 
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Thus, for example, the sixth time derivative of the output of the nonlinear block 2 , namely has the value of 
30.31 at t = 1. The algorithm is shown as a block diagram in figure 2.2. Clearly the process could be continued 
if there were additional polynomial nonlinearities in the sequence, so multiple nesting of polynomials is easily 
handled by means of dynamic forms. 


a = 


( 1 1 -n 

11 — - — 

2 3! 4! 


b = (l -0.5 -0.5 O.l) 



Figure 2.2. Scalar forms for the example of nested polynomials. 


In the algorithm poly' n the coefficients are fixed. We denote polynomials with time-variable coef- 
ficients by poly n . Thus 

z = poly n (a, x) = a.Q + a\x + . . . + a n x n ( 2 - 22 ) 

is a polynomial with possibly variable coefficients. The corresponding algorithm for dynamic forms is 
obtained from the poly' n (a, x) algorithm in equation (2-20) by replacing not only z and x by their forms 
but also aj by SF m (ai ) for i = 0, . . . , n: 

algorithm : SF m (z ) = poly n [S F m (a) t S F m (x)] 

SF m (z) = SF m {a n )*SF m (x) 

do i — 1 , n — 1 

SF m (z) = [ SF m (a n _i ) + SF m (z )] * SF m {x ) 

end do 

SF m (z) = SF m (a 0 ) + SF m (z) 

Nonlinear functions of several variables are frequently given in practice in tabular form. Thus, 
for example, the subroutine generating the aircraft total force and moment may contain a table and 
interpolating routine representing the functional dependence of the moment coefficient on the angle of 
attack and the Mach number, C m = f( a > Mach). Sometimes such tables may be represented by nested 
polynomials, which for two variables may take the following form: 

2 = f(x, y) = (&oo + • • • + &o n y n ) + (&io + • • • + hny n )x + • ■ • + (ho + • • • + h n y n ) xk 
The algorithm for computing m time derivatives of / is easily constructed: 
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algorithm: SF m (z) = f[SF m {x), SF m (y)} 

do i = 0, k 

SF m (a i )=poly n [bi, SF m (y)} 

end do 

SF m (z) =poly k [SF m (a),SF m (x)] 
where bi = ( fyo. • • • . b{ n ). 

Thus we are beginning to compute easily time derivatives of fairly complicated functions. In the 
preceding example an algorithm was constructed from simpler algorithms that were built up from still 
simpler functions. At the bottom of such a hierarchy are functions from a basic set that already includes 
sum and product. Next, more functions are added to this set. 


Inverse, Division, Fractional Powers, Logarithms, and Exponentials 


Many standard functions may be converted to functions of scalar forms by means of the Leibnitz 
product rule. In this section we derive the algorithms for x _1 , x/y, x°/*\ e x , and lnx. Trigonometric 
functions cosx, sinx, and arctan(y,x) will be derived in the next section. For convenience, we shall 
denote, say, SF k [f{x, y, . . .)] by f[SF k (x ), SF k {y ), . . .]. Thus, for example, 

ln5F m (x) = SF^lnx) = (lnx, (lnx)* 1 ), . . . , (lnx)^ ) (2-23) 

The notation will allow such easy visual and machine translation as, for example, 

z = g— ax 2 — 6y 2 ^ SF m ^ = e -a[SF™(z)] 2 -6[SF™(y)] 2 

In the following discussion it is understood that for any x, /[SF m (x)] is defined only if /(x) is defined. 


If z = x/y and j/ / 0, then, since x = zy, the time derivatives are given iteratively by the 
following algorithm: 


z(°) = xy 1 

< 

z (k) = ( ' x (k ) _ £fc =i Qk z { k -i)y(i)^ y - 1 if 0 < k < m 

This division algorithm will be denoted as 

SF m {z ) = SF m (x)/SF m (y) 

It may be noted that the cancellation law holds for forms: 

[SF m (w) * SF m (x)]/[SF m (w) * SF m (y )] = S F m (x) / S F m (y) 

In particular, if z = y _1 , then, since 2 = 1/y, the division algorithm produces the inverse. 

z(°) = y -1 

< 

z ( k ) = - (£* =1 cf 2 ( fc - f )y(0) if 0 < k < m 


(2-24) 


(2-25) 

(2-26) 


(2-27) 
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The effect of this algorithm will be denoted as 

SF m (z) = [SF 771 ^)]" 1 (2-28) 

and referred to as the scalar form inverse. Note that 

SF m (x ) * [SF m (x)}- 1 = [5F m (x)] _1 * SF m (x) = SF m ( 1) 

If z = x a ! b for constant integers a, b with b > 0 and x ^ 0, then bxz = azx\ hence, on application of 
the Leibnitz rule to both sides, 

k k 

£ C*bx®z( k - i > = CfazWxl*-*) (2-29) 

i=0 i~0 

But since z^ k ~^ = z^ k ~ ' i+1 ), and similarly for x, 

k k 

Ci bx^z^ k ~ i+l ^ = Cf az (7 M fc_i+1) (2-30) 

?=0 i— 0 

Or, peeling off the first (i = 0) term, for 1 < k, 

bxz^ k+V > -(- £ CfbxWz( k ~ i+ V = az(°M fc+1 ) -|- £ Cfaz(*V*- i+1 ) (2-31) 

i=l i=l 

Consequently, we obtain the following basic algorithm 
r z(°) = x a > b 

< z = (bx) - ^^ 1 ) (2-32) 

k -z(k +1 ) = (6x) _1 {az^x^ +1 ) -I- Yli=i C k [az^x( k ~ i+1 } — &x(*)z(* — * +1 )]} for 1 < k 

The effect of this algorithm will be denoted as 

SF m (z) = [SF m (x)] a / b (2-33) 

We note that 

[SF m (x)) a / b * [SF m (x)] c / d = [SF m (x)}( a M+W (2-34) 

If z = x”, with i^0 and n > 0, then 

SF m (z) = [SF^x)] 71 / 1 (2-35) 

which will be simplified to 

SF m (z) = [5F m (x)] n (2-36) 
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It may be noted that this algorithm involves only one Leibnitz convolution instead of n as in the simple 
algorithm discused previously. 

If z = |x| then, of course, 

SF m {x) if x > 0 

|SF m (x)| = J (2-37) 

if x < 0 

If z = e x , then z = zx and so, for 0 < k < m — 1, the Leibnitz product rule gives 

z (k+\) = £ cfzP-VxQ+v 

i=0 

We denote this algorithm as 

SF m (z) = e SFm W 


x 1 x^ 1 ) and for 


(2-40) 


(2-41) 
(2-42) 

Other functions may now be built up. For example, the scalar form of a Gaussian, z = be x , is 
computed easily by first calling the dynamic-form squarer, then the sign change, and then the exponential: 

SF m (z) = 6e-[ 5Fm ( x >] 2 (2-43) 

Note how the notation keeps track of the calling sequence. In effect the result of each call has an 
independent existence. On the other hand, the algorithm SF m (z) = SF m (be~ x ), while true overall, 
is not parsed as a calling sequence of independent subroutines. 


Conversely, if z — lnx for x > 0, then z^ = lnx and since xi = x, z = 
l<fc<m— lthe product formula leads to 


*(*+!) = X ” 1 x( fc+1 ) - £ cfxW*(*- i+1 ) 

V i= 1 


This algorithm will be denoted as 


It may be noted that 


SF m (z) = \nSF m {x) 


e In[SF™(x)] = SF m, x j 


(2-38) 

(2-39) 


Trigonometric Functions 

To obtain trigonometric functions, it is convenient to group cos and sin as follows: 



(2-44) 
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Then, since 


i = xQz 


(2-45) 


where 



the derivatives of cos x and sin x are given by 


z (k+l) = £ 

i=0 

The algorithm (2-47) for the scalar form of the cos and sin functions will be denoted as 




^ / cos SF m (x)' 

\ sin SF m (x) 


(2-46) 


(2-47) 


(2-48) 


Conversely, suppose that z is redefined as 


z = r 


( cosx 
sinx 


(2-49) 


with r > 0 so that 


x = arctan(22, z i) 


(2-50) 


Then the derivatives of x may be computed without assuming that the norm of z is 1, as follows: From 
equation (2-49), r 2 = z T z and, upon differentiating equation (2-49), i = rr~ l z + xQz. Premultiplying 
the last equation by z T Q and noting that z T Qz = 0 and QQ = —I, we obtain 


rp 

xz z = z Qz 


(2-51) 


Let 


Then 


Now we raise the order. 


U = Z T Z = Z\Z\ + Z%Z2 
i 

L w = z T Qz = -21^2 + %2 Z 1 


X = w/u 


(2-52) 

(2-53) 


and 


f SF m (u) = SF m (zi) * SF m (z i) + SF m (z 2 ) * SF m (z 2 ) 

\ SF m - 1 (tn) = -SF m-1 (ii) * SF m - l (z 2 ) + S , F rn - 1 (i 2 ) * SF m ~ l (z\) 


(2-54) 


SF m ~ l {x) = SF m ~ 1 (w) /SF m ~ 1 (u) 


(2-55) 
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Equations (2-50, 2-54, and 2-55) define the arctangent algorithm, which will be denoted as 

SF m (x) = arclan[SF m {z 2 ),SF m (z 1 )] (2-56) 

Of course, if it is known that r = 1 identically, that is SF m (u ) = (1,0), then 

= arctan(z2, zi) 

< 

k = w^\ 0 < k < m — 1 


If a nonlinearity is modeled by a few terms of its Fourier series (here k is the fundamental wav~ number), 

n 

z = fourier n (k, a, b,x) = a o + X! a i cos (ikx) + sin (ikx) (2-57) 

i=l 

then SF m (z ) is given by 


fourier n [k, a, b, SF m (x )] = aoSF m (l) + a i cos[ifc5F m (x)] + 6j sin[tA:5F m (x)] 


i=l 


(2-58) 


Functions of several variables can be treated similarly. Multivariable polynomials find use in 
modeling aircraft force and moment generators (ref. 24). Multidimensional Fourier series and other 
functions may be used for describing winds, terrain, or other fields that influence the aircraft. For 
example, if the generating function is a polynomial in p variables each of maximum degree n, 

z= ^2 a(i \ , . . . , i p )x i 1 . . . Xp p (2-59) 

then SF m (z ) may be evaluated (assuming all Xj > 0) by the following algorithm: 


do i — l,p 

SF m ( yi ) = In SF m ( Xi ) 

end do 


SF” l (z) 







The following example illustrates the application of dynamic forms. 
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Example 2.2. As an example of the application of scalar forms, consider a nap-of-the-Earth (NOE) problem 
(ref. 25). A simple version of this problem may be stated as follows: Suppose that a plan view of the trajectory to 
be flown by the helicopter is given. The altitude must be maintained at a fixed, usually minimal, distance above 
the terrain. The problem is to determine whether the resulting trajectory is within the helicopter performance 
limits. The intermediate problem is to determine altitude and its time derivatives. Suppose that the horizontal 
velocity of a helicopter is given in polar coordinates by horizontal speed s(t) and heading ip(i). Then the plan 
view of the path evolves according to 

n-n 

\y / \ s sin ip / 

and , 

/x\ /i(0) + / 0 i\ 

V y J V y{ 0) + Jo y ) 

where the initial condition x(0), y{ 0) is also given. Suppose, in addition, that there are some irregular vertical 
obstructions that are collectively covered by a Gaussian cap: 


h — hr 


-( x/tJ x Y-{y/<r y y 


where h is the altitude. The helicopter is required to stay on that cap while flying the given horizontal trajectory. 
The problem is to determine the vertical speed and the next four of its derivatives. It is easily done; simply 
replace the variables by their forms, as follows: 


5F 4 (x) 


SF 4 (y) 


SF 5 (x) 

SF 5 (y) 


/SF 4 ^)* cos SF 4 (ip)\ 

\ SF 4 (s) * sin SF 4 (ip) / 
x(0) + Jg SF°(x) \ /SF*(x) 

y(0) + tiSF°(y))\SF 4 M 


SF*(h) = 

The scalar form SF 5 (h) contains the altitude h and five of its time derivatives. These higher derivatives of 
h are needed for testing constraint satisfaction, such as acceleration h, which affects pitch angle and power 
requirements, h,( 4 \ which has a direct bearing on moment requirements, and h^ 5 \ which affects the usually 
limited control rates. Details of the force and moment model would be needed to compute the actual values of 
controls. Of course it is possible to compute other functions of scalar forms. For example, the total translational 
kinetic and potential energy divided by weight is 

e =-^r +h 

where g is the acceleration of gravity. The algorithm for computing translational power and power rate is obtained 
by simply replacing variables by their forms: 


SF 4 (e) _ + SF \h) 


A numerical example is given in appendix A. 
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Partial Derivatives 


We have developed in the preceding sections several functions of dynamic forms. With these 
functions as building blocks many other practically useful functions can be assembled. For example, 
suppose that there is an algorithm for computing 

z — f{x , u ) (2-60) 

for scalar z, x , u, and that this zero-order algorithm has been converted to order m: 

SF m {z) = f[SF m (x),SF m {u )] (2-61) 

as shown in figure 2.3. The input consists of the scalar forms SF m (x ) and SF m (u), and the output 

is the scalar form SF m {z ) = f[SF m (x), SF m (u)}. Now we proceed to compute, using only this 

algorithm, several useful objects associated with /. 

Consider the partial derivatives of /. Since 2 = f(x,u). 


z = f x x + f u u 


(2-62) 


where f x denotes the partial derivative of / with respect to x. Choose an input to be 


' SF m (x) = (x,l) 
^SF m (u) = (u, 0) 


(2-63) 


where zero padding is assumed: ( x, 1 ) = ( x, 1, 0, . . . , 0 ) and ( u, 0 ) = ( u, , 0, . . . , 0 ). Then the value 
of the re-partial derivative of / at x is given by i, that is, by the content in the number one location of 
SF m (z) 

fx{x, u) = z = {5F m (2)}i = {/[(*, 1), (u, 0)]}i (2-64) 

where we use the notation 

f if 0 < k < m 

{SF m (y)} k = 

{ 0 else 

The purpose of equation (2-64) is not to replace a familiar and convenient notation by one that is obscure 
and awkward; the purpose is to show that the useful object, the rr-partial derivative of / at (x,u), may 
be computed by means of scalar forms as the value in location one of the output of the scalar form 
algorithm / whose input is (x, 1) and (u, 0). 


X ^ 

f 

SF m (x) ► 

► z => 

f 

u — ► 


SF m (u) ► 



Figure 2.3. Algorithm for a composite function /. 


17 



(2-65) 


The value of the k th partial of / at x with respect to x is given by 

f x...x = 1), (li, 0)]}*, 1 < k < m 

Similarly, 

fu...u = {/[(x,0),(u, 1)]}*, 1 < fc < m 


( 2 - 66 ) 


Furthermore, since 


/ 


fxxii + fxX 2fxuXU + /itii -f /uuirM 


and 

fxx = {/[(x,l),(u,0)]} 2 

fuu — {f[(x, 0), (u, 1)]} 2 

it follows that 

fxu = f[(x , 1), («, 1)] - f[(x, 1), (u,0)] - /[(x,0), (u, 1)]} 2 


(2-67) 

( 2 - 68 ) 


(2-69) 


Function Inverse 

The algorithm for computing the dynamic form of the inverse of a function is often of practical 
interest. Suppose that an algorithm is given for computing a possibly time-varying function 

z = f(x, u , t) (2-70) 

for scalar z, x, u and time t; that this zero-order algorithm has been raised to order m: 


SF m (z ) = f[SF m (x), SF m (u), SF m {t )] (2-71) 

and that we have the (partial) inverse / -1 of / 

u = f~ l (x,z,t) (2-72) 

so that for all x, u, z, t 

f[x,f~ l (x,z,t),t] = z (2-73) 

Often in practice / -1 is obtained numerically by an algorithm such as Newton-Raphson. The extension 
of f~ l (x, z,t) to 

SF m (u) = r 1 [SF m (x),SF m (z),SF m (t)] (2-74) 

producing not only u but also m of its time derivatives, may be computed as follows: The time derivative 
of equation (2-70) is 

Z = fxX + fuU + ft (2-75) 

But, according to equation (2-71), 

i = {/[(x, x), (u, u), ( t , l)]}i = {f[(x, x), (u, 0), (<, l)]}i + f u u (2-76) 
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Hence it follows that 


U=(fu) 1 (z - {/[(*, x)>( u .°MM)]}l) (2-77) 

where the inverse of the partial derivative f u is 

(/«)“* = {/[(l,0),(u,l),(t,0)]}r 1 (2-78) 

Similarly, 

Z = {/[(*, X, x), ( u , u, 0), (t, 1 , 0)]>2 + fuU (2-79) 

so that 

U = (fu)~ l (z - {/[(x, X, x), (u, u, 0), (t, 1, 0)]>2) (2-80) 

and for higher derivatives, 

= (fu)- 1 (z (fe) - {f[(x, . . . ,x^ k )), (u, . . . ,^^,0), (t,l)]} k ) (2-81) 

Thus, in addition to the base point u = / _1 (x, z), only derivatives of / are needed for the construction 
of derivatives of / -1 . Hence the algorithm: 

ALGORITHM: SF m {u) = f~ l [SF m (x), SF m (z)} 

SF m (u) = (f~ l (x,z, *),0 0) 

(/«) _1 -{/[(*■ 0),(«,l),(t,0)]}f l 

do k = 1,771 

« (l) = ( A )" 1 (z (k) - {f[SF k (x), SF k (u), SF*(()]} t ) 

end do 

That is, first SF m (u ) is loaded with the base point u = f~ l (x,z,t ) and zero derivatives; then the 
inverse of f u is computed as in equation (2-78); finally, equation (2-81) is iterated from 1 to m. The 
combined action of / -1 and / is shown in figure 2.4. The input is SF m (x), the desired evolution is 
SF m (z), and the required control is SF m (u). So far we have been considering static systems. Next 
we consider dynamic systems such as differential equations. 


SF m (x) 

SF m (z) 


I 


t —1 


SF m (u) 


1 


SF m (z) 


Figure 2.4. Useful factors of identity. 
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Differential Equations 

Suppose that we are given a possibly time-varying, first-order scalar system 

x = f(x,u,t) (2-82) 

where x is the state and u the control, and we wish to compute derivatives of x up to order m given 
the derivatives of u. First, as described previously, we convert the zero-order algorithm 

Z = f(x,U,t ) 


to order m: 

SF m (z) = f[SF m (x), SF m (u), ( t , 1)] (2-83) 

where, as before, SF m (t ) = (£, 1, 0, , 0) is abbreviated as (£, 1). Then 

*<* +1 > = {/[SF*(z), SF k (u), (t, 1)]}* (2-84) 

Consequently the time derivatives of x can be computed iteratively from the initial condition x and the 
control S'F m (w): 

ALGORITHM: [5F m (a:)] = Tf[ X, SF m (u )] 

SF m (x) = (x, 0, . . . , 0) 

do k = 0, m 

*<*+!> = {/[SF*(:r),SF*( U ),(t,l)]}* 

end do 

That is, first SF m {x) is loaded with the base point x and zero derivatives; then the derivatives x^) are 
obtained iteratively by means of equation (2-84). The effect of the algorithm will be denoted as 

[SF m (x)] = T f [x, SF m (u)] (2-85) 

where the subscript in Tf is a reminder that the function / raised to order m must be provided. 

Now that we have SF m (x), we can approximate the solution <j>(x,t + r) of equation (2-82) by the 
first m + 1 terms of its Taylor series: 


TYl 

4>{x, t + r) « tay[SF m (x), r] = ^ x^r k /k\ 

k=0 


( 2 - 86 ) 


The Tf algorithm may be easily generalized to scalar differential equations of higher order. Thus, 
suppose that in equation (2-82) x € R n , u € R, and that we have zero-order algorithms 

Zi = fi(xi,...,x n ,u,t) (2-87) 

for 1 < i < n and these algorithms have been converted as discussed previously to 

SF m ( Zi ) = fi[SF m (x i), . . . , SF m (x n ), SF m (u), (£, 1)] (2-88) 
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(fc+1) ( k ) 

Then, since ' = z\ ’ for i = 1, . . , , n, the time derivatives of the n-dimensional state x may be 
constructed from the base point x and the control SF m (u) as follows: 

algorithm: [5F m (a:)] = T f[x , S F m ~ l (u)} 
do i = l,n 

SF m (xi) = (xi, 0 ,..., 0 ) 
end do 
do k — 0, m 
do i = 1,71 

zf +1) = (fi\SF m ( X1 ) SJ^(*„),SF™(u),(f,l)]}* 

end do 
end do 


The input is the state x and the evolution of control SF m (u); the output is the evolution of the 
n state coordinates SF m (x{). In the algorithm the base point x with zero derivatives is loaded into 
SF m (x ); then the derivatives are computed by columns using /j. It may be noted that, in general, the 
control u must be available to order m — 1 for the computation of Xj to order m. An important special 
case occurs for pure feedback systems for which the function / above has the following triangular 
structure: 

fi(xi,...,x n ,u,t) = fi(xi,...,x i+ i,t), 1 <i<n (2-89) 

so that the state equation has the following form: 


( x\ \ 

X2 

X3 


fl(xi,X 2 ,t ) \ 

f 2 (xi,x 2 ,x 3 ,t) 
fs(xi, x 2 , X3, X4, f) 


(2-90) 


\Xn/ \f n (x 1 ,X2,...,x n ,u,t)/ 

For such systems u to order m—n is sufficient to compute X{ to order m+l — i. In particular, derivatives 
of it are not used for the computation of x\ to order n, and u itself is not used in the computation of 
x\ to order n — 1. The computations in Ty for n = m = 4 may be represented as follows: 


- 

JO) 

x l 

JD 

x l 

x (2) 

x \ 

J3) 

x l 

J4) 

x \ 

J 5) 

x l 

J6) 

. f . 

X\ 

• 

* 

* 

* 

* 

O 

0 

X2 

• 

* 

* 

* 

0 

O 

0 

^3 

• 

* 

* 

0 

0 

O 

0 

X4 

• 

* 

O 

0 

0 

O 

0 

u 

• 

O 

O 

0 

0 

O 

0 


where • denotes the initial data, and * and o denote significant and insignificant new entries, repectively. 
The computation flows from left to right starting with the supplied first column. The algorithm Tj 
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specialized to the case in which the output of the system is x\ will be denoted as T\ j, that is, 

SF m (x i) = T lf l(x,SF m - n (u)] (2-91) 

where, if m < n, then the entries in the scalar form of u are not used by the algorithm. The block 
diagram of 7\ y, which we will refer to as the forward solution since u(t) is the input and the Taylor 
series of xi(t) is the output, is shown in figure 2.5. The input consists of the base point x and the scalar 
form SF m ~ n (u ) of control u. The output is the scalar form SF m (x i) of the lowest state coordinate 
xi. Only the first m — n derivatives of u are used by the algorithm. If m<n, then u is not needed at 
all in the computation of x^ m \ This condition, in which no entries of a scalar form SF m (x) are used 
by the algorithm, will be denoted in diagrams by SF m (x ) = 0. 


x 


SF m_n (u) 



SF m ( Xl ) 


Figure 2.5. Forward solution. 


Example 2.3. Suppose that n = 4, the state equation is 

X2 \ 

a sin(6xi)x2 + [2 + cos(c<)]x3 


/i 


X2 




\X 4 / 



X4 


\ 


U 


I 


and a, b, c are constants. The first step is to raise fi to order m as follows. 

A = SF m (x 2 ) 

h = osin[65F m (xi)] * SF m (x 2 ) + {2 + cos[cSF m (t )]} * SF m (x 3 ) 
h = SF m (x 4 ) 
h - SF m (u) 

Then the forward algorithm T\ f is constructed. 


(2-92) 


(2-93) 
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ALGORITHM: SF m (x ) = T X f[x, SF m " 4 (u)] 
do 2 = 1,4 

SF m ( Xi ) = (xi, 0 — , 0) 

end do 

do k = 0, m — 1 

x{ k+l) = {SF"-(x 2 )H 

r- 2 k+1> = {osin[6SF m (x 1 )] * SF m (x 2 ) + {2 + cos(cSF ra (i)]} * SF m (x: s )‘ lk 
4* +l) = {SF™(x t )) k 
x<‘ +1) = ’,SF™(u,} k 

end do 


This very simple example, for which even hand computation of derivatives is relatively easy, is 
intended only to illustrate the technique. 


Dynamic Inverse 

Consider again the first-order, scalar, dynamical system in equation (2-82). It is useful to be able 
to compute what the control signal must be so that the system output will track a given function of time. 
The evolution of control u that will produce the desired evolution of the output x may be obtained as 
follows: Construct the inverse u = ) of z = f(x,u,t) so that 

f[xj~ l {x, Z,t),t] = 2 

and raise its order to m by means of the function inverse algorithm as described previously. 

5F Tn (tx) = r 1 [5'F TO (x),5F m ( 2 ),(f,l)] (2-94) 

Then the control evolution SF m (u ) that will produce the desired output SF m (x ) is obtained by imposing 
the constraint z = x as given by the following algorithm: 

ALGORITHM: [x, SF m ~ n (u )] = Tjy 1 [SF m (x)] 

SF m ~ 1 (z) = SF m ~ l (x) 

SF m ~ l (u) = f~ 1 [SF m ~ 1 (x), SF m ~ 1 (z), (t, 1)] 
x = {SF m (x)} 0 

That is, first SF m (x ) is shifted to the right and loaded into SF m ~ l (z ) to get 

( 2 ,...,^ m - 1 )) = (x( 1 ),...,x( m )) 

Then a call to the f~ l algorithm produces SF 771-1 (w), and the state x is just the zero term of SF m (x). 

This algorithm may be generalized to higher order differential equations, provided that they are of 
the pure feedback form as in equation (2-89), and that, for 1 < i < n, each fj is invertible with respect 
to Xj + i and f n is invertible with respect to u (we denote these inverses by /“*), 

Zi+l = • ,Zi,t) for 1 < i < n 

< (2-95) 

k u = /j _1 (a;i,..., 2 t+l,f) fori = n 
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Thus, 


ALGORITHM: [x, SF m n (u)] = T^ 1 [SF m (x \ )] 
do k = 1, n — 1 
SF m ~ k (zj c ) = SF m ~ k (x k ) 

SF m ~ k (x k+ 1 ) = f^lSF^-^xi) SF m ~ k (x k ), SF m ~ k (z k ),(t, 1 )] 

end do 

SF m - n (z n+l ) = SF m - n (x n ) 

SF m ~ n (u) = , SF m - n (x n ), SF m ~ n (z n+1 ), (t, 1 )] 

do i = 1, n 

Xj = {SJ™(*j)}o 

end do 


Thus, the inversion is accomplished downward by rows until all n coordinates of x have been 
computed; then a call to /„ 1 produces the control and its derivatives; finally the base point x is 
assembled from the zero-order terms. The state must be at this value of x in order for the evolution 
SF m (x i) to be possible. The flow in algorithm T~ f l for n = m = 4 may be represented as follows: 


I 

JO) 

X 1 

x (1) 

Xl 

x (2) 

Xl 

x (3) 

Xl 

X {4) 

X 1 

x (5) 

Xl 

x {6) 

Xl 

Xl 

• 

• 

• 

• 

• 

o 

o 

*2 

* 

* 

* 

* 

0 

o 

o 

X3 

* 

* 

* 

o 

0 

o 

o 

X4 

* 

* 

o 

o 

o 

o 

o 

u 

* 

o 

o 

o 

0 

o 

o 


The computation proceeds downward from the initially supplied first row. In the usual application of 
linearization techniques to the control of pure feedback systems, m = n so that only the control u is 
obtained. Cases with m > n are of interest when derivatives of u are needed for path planning. The 
block diagram of this algorithm, which we will call the inverse solution since the input is x\ (t) and 
the output is u(t), is shown in figure 2.6. The input is the scalar form SF m (x\) of the lowest state 
coordinate xi. The output consists of the complete state base point x = (xi,X 2 , . . . ,x n ) T and the 
scalar form SF m ~ n of the control variable u. 


SF m (x n ) 



x 

SF m_n (u) 


Figure 2.6. Inverse solution. 


% 

If m <n, then neither u nor any of its derivatives are computed. 
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Example 2.4. The inverse f i 1 of in example 2.3 is easily computed: 

/f 1 = SF m ( Zl ) 

/- 1 = { SF m (z 2 ) - asin[65F m (xi)] * SF m (x 2 )}/{2 + cos[cSF m (<)]} 
f -l = S F m (z 3 ) 


(2-96) 


/ 4 _1 = SF™(z 4 ) 


Equation (2-96) defines the inverse algorithm 7^ j. 

ALGORITHM: [x, 5J rm_4 (u)] = [SF m (xi)] 

SF m ~ 1 (z 1 ) = SF™- 1 ^!) 

SF rn - 1 (x 2 ) = SF m - 1 (z 1 ) 

SF m ~ 2 (z 2 ) = 5F m-1 (x2) 

SF m ~ 2 (x 3 ) = {SF m-2 ( 2 2 ) - a sin[6SF m-2 (xi)] * SF m ~ 2 (x 2 )}/ {2 + cos[cSF m - 2 (t)]} 
SF Tn ~ 1 {z 3 ) = SF m ~ l (x 3 ) 

5F m ~ 3 (x 4 ) = SF rn ~ 3 (z 3 ) 

SF m ~ 1 (z 4 ) = SF m ~ 1 (x 4 ) 

SF m ~ 4 (u ) = SF m - 4 (z 4 ) 
do i = 1,4 
Xj = {5F m (xj)}o 

end do 


The example illustrates how easy it is by means of dynamic forms to propagate the derivatives 
backwards through the system. 

It may be noted that T^} is the inverse of T\ j in the following sense. For any evolution SF m (x 1 ) 
of the output variable x\, 

= SF™(x 1 ) (2-97) 

This significant relation is shown as a block diagram in figure 2.7 An application of this relation to 
automatic control is obtained by inserting the plant between and T ) j as discussed next. 


SF m ( Xl ) 



SF m ( Xl ) 


Figure 2.7. Useful factors of identity. 
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Automatic Control 


Consider the problem of tracking a given reference input in the presence of disturbances. Let the 
plant be described as follows: The state x p € R n , the control u p € R, the output y = x\ p , and the state 
equation is possibly nonlinear and time-varying but of maximal relative degree (pure feedback form): 

/ fl ( x lpi x 2pi t) \ 


Xp — f ( X p , Up, t ) 


/ 2 (:rip,X2p,X3p,t) 


(2-98) 


V fn (% lpi x 2 p> • ■ ■ i %np-> u p, t) ) 


The structure of a possible control-system design (ref. 26) is shown in figure 2.8. There are 
three subsystems in addition to the plant, namely transformer, regulator, and command generator. The 
command generator provides at every t not only the required value of the commanded output x\ c but 
also at least n of its derivatives, i.e., SF n (x i c ). The plant state x p is passed through 7 \ j to obtain 
SF n ~ l (xi p ). Note that this step does not involve any differentiation of signals that are corrupted by 
noise. Furthermore, no element of SF(u) is needed in this step. Next the error form is computed to 
order n — 1 in the regulator by comparing the output form with the input form: 

SF n -\e) = SF n ~ l (x ip) - SF n ~ l (x lc ) (2-99) 



GENERATOR REGULATOR TRANSFORM PLANT 

Figure 2.8. A structure for asymptotic trackers with pure feedback plant dynamics. 
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Then the order of the error form is raised by one by the regulator law k, and the error form is added to 
the reference, thereby raising the order by one of the desired output form: 

SF n (x ip) = SF n (xi c ) + SF n (e) (2-100) 

Finally, the desired evolution of the output is passed through the inverter T^ 1 , whose output SF®(u) 
drives the control u p . 


It may be noted that the transformations T^ 1 and 7\ j make the plant look like an (invariant) string 
of integrators, n integrators long. Thus, the plant has been transformed into a Brunowsky canonical 

form (ref. 14) with the Kronecker index equal to n, where (x \ p , . . . , x'^ 1 ) is the transformed state 
(n) 

and x\ is the transformed control. Consequently, the system may be regulated by means of a simple 
linear law 


e 


n — 1 


i = 0 


( 2 - 101 ) 


despite changes in the perturbation model of /. In effect, the transformations 7 \ j and T^jr provide 
automatically the necessary gain scheduling and therefore the system output x\ p will track asymptotically 
any input x\ c that is differentiable n times. The derivatives of the input must be available, and the 
plant dynamics must be of the pure feedback type. It should be noted, however, that only the state is 
needed and none of its derivatives. Noise is not an issue. The control algorithm denoted as Cf may be 
summarized as follows: It is assumed that Cf is called at the implementation sampling rate, which is 
sufficiently fast relative to the dynamics of the plant. 


ALGORITHM: U p = C f[x p , 5i ?n (xi c )] 

SF n -\x lp ) = T lf (x p ,<&) 

SF n ~ l (e) = SF n ~ l (x lp ) - SF'-Hx ic) 

e<") = Jfc[SF n “ 1 (e)] 

SF n (x lp ) = SF n (x ic) + SF n (e) 
[5F°(wp),x p ] = T^f[SF n (xi p )] 
u p = {SF°(u p )} o 


The input to the algorithm consists of the (estimated) plant state x p and the generator command 
SF n (xi c ). The output is the plant control u p . In the algorithm the symbol 0 indicates that the content 
of SF n ~ l (u) does not matter for this call of the subroutine T\ f. The output x p (next to last line in the 
algorithm) is the plant state. 


It may be noted that Cf is a general algorithm in the sense that it applies to any plant defined by 
/ such that 1) / is of the pure feedback type and 2) / is a composite of the functions that have been 
elevated to order m. In such a case the zero-order function / can be automatically raised to order m, at 
which point all the steps in Cf become fixed. A numerical example illustrating the operation of control 
systems with the structure shown in figure 2.8 is given in appendix B. 

In summary, we have shown that scalar forms may be used to organize effectively the design of 
automatic control systems for the class of plants having pure feedback dynamics /. The key step is to 
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raise the order of / from zero to n, but that step is made routine by the formalism of dynamic forms 
in which elementary functions and operations have been translated into their dynamic -form equivalents. 
In the remainder of the report we translate additional functions frequently occurring in the design of 
control systems. Further discussion of control algorithms, including plant dynamics with transmission 
zeros, will be presented in a future report. 

3 VECTOR AND MATRIX FORMS 

We are concerned mainly with three-dimensional vectors. Vectors and dot and cross products will 
be denoted as usual. Thus, x and y are vectors, and x ■ y and x x y are their dot and cross products. 
Right-handed orthonormal axis systems will be denoted by a double arrow. Thus *a is an axis system 
( a\ a 2 < 13 ) where, for i,j = 1,2, 3, a; • aj = 6ij and 03 = a\ x 02 - The a-coordinates of x will be 
denoted by the column matrix x a . Thus, 

3 

X ^ ] ^aj^j ( 3 - 1 ) 

3 = 1 

the dot product x ■ y = x^y a , and the cross product z = x x y 

z a = S{x a )y a (3-2) 

where, for any x £ /? 3 , the skew symmetric matrix 

0 x 3 -x 2 
S(x) = —£3 0 xi 

V x 2 — x\ 0 

The transpose is denoted by (-) T and <5j denotes a column with 1 in row i and zeros elsewhere. 

Next, consider vector forms. Let *a be an axis system. The vector form of a vector v in *a is 
defined by 

VF" l (v a ) = ( Va , ia 4“)) (3-4) 

where v a £ R n . The zero and unit forms will be denoted by 

VF m ( 0 ) = (0 ... 0) (3-5) 

and 

VF m (6 l ) = (6 i ... 0) (3-6) 

where 5j is a column of zeros except in row i, which contains 1. We denote the scalar form of the i th 
component of a vector form V F m (v a ) with a subscript as follows: 

SF™(v ai ) = Vf-Mi = {v ah v ai v<?) (3-7) 
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Next consider matrices. A matrix form of order m is defined as a matrix and rn of its time 
derivatives. 


MF m (A) = (a,A,.., 

with A G R n x R n . Zero and identity are given by 


(3-8) 

MF m { 0) = (0 0 . 

... 0) 

(3-9) 

MF m (I) = (/ 0 , 

... 0) 

(3-10) 

and the transpose 

[MF m (Z)] T = MF m (Z T ) 

(3-11) 

The scalar form corresponding to the Zjj element of Z will be denoted as 


SF m (z i j) = MF m 

{Z)ij 

(3-12) 


Algebra 

Several useful functions of vector forms and matrix forms follow easily either from already-defined 
functions of scalar forms or from the Leibnitz rule. 

If*a = ax a + by a for either scalars or matrices a, b constant in time, then, of course, 

VF m (z a ) = aVF m (x a ) + bVF m {y a ) = VF m (ax a ) + VF m (by a ) (3-13) 

If z a = ax a and a is a scalar function of time, then the derivatives of z a are given by 

zi k) = tc^ k -^ (3-14) 

i=0 

We denote the resulting algorithm as 

VF m (z a ) = SF m {a ) * VF m {x a ) (3-15) 


If z = x ■ y, then 


which will be denoted as 


If z = x x y, then 


■S k) = E 


i = 0 


SF m (z) = VF m {x a ) • VF m (y a ) 


z { a k) = -J2CiS(x { a i] )y^ 


i = 0 


(3-16) 

(3-17) 


(3-18) 
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which will be denoted as 


VF m (z a ) = VF m (x a ) x VF m (y a ) 


( 3 - 19 ) 


If Z = aX for time varying scalar a and matrix X, then by the Leibnitz rule 

k 

z i k ) ~ Cfa^xW 

i=0 

So we have an algorithm for multiplication of a matrix form by a scalar form. 


MF m (Z) = SF m (a ) * MF m (X) 


If z a = Xy a with matrix X and vector y a , both time-varying, then 

4 k) = £ CfxMyP 

i=0 

which is an algorithm for time-varying transformations to be denoted as 


VF m (z a ) = MF m (X ) * VF m (y a ) 


If Z — XY with time-dependent matrices X and Y, then, again by the Leibnitz rule, 


z (k) _ J2 CiX^-^Y® 

i = 0 


This algorithm will be denoted as 


MF m (Z ) = MF m (X) * MF m (Y) 


If X is invertible, then in direct analogy to scalar forms the derivatives of X 1 are 
Z(°) = X -1 , and for 1 < k < m 

Z^ = - ^ cf Z^-^X^ j z 

So we give meaning to 

MF m (Z) = [MF m (X)}~ 1 

It is now possible to construct higher level functions of dynamic forms. For example, if x 
[ujua] 1 / 2 , then the evolution of the magnitude of v a is given by the algorithm 

SF m (x ) = \VF m (v a )\ = [ VF m (v a ) • VF m {y a )} 1 ! 2 
where the operations on the right side are all defined. 

The unit vector form corresponding to u parallel to v in a* is given by 

VF m (u a ) = VF m (va)/\VF m (v a )\ 


( 3 - 20 ) 

( 3 - 21 ) 

( 3 - 22 ) 

( 3 - 23 ) 

( 3 - 24 ) 

( 3 - 25 ) 
given by 

( 3 - 26 ) 

( 3 - 27 ) 

= H = 

( 3 - 28 ) 

( 3 - 29 ) 
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Derivatives 


Suppose that x,z e R n and 2 = f(x) and that we have raised the order of the vector function / 
from zero to m, 

VF m (z) = f[VF m {x)] (3-30) 

Consider constant vectors x and y such that y T y = 1. Let 

z = f(x + ty) 

Then the derivative of / at x in the direction of y at t = 0 is 

2 = fx(x)y 

So by choosing 

VF m (x) = (x,y,0,...,0) = (x,y ) 
we obtain the Taylor series expansion of / at x in direction y. 

vn 

f(x + sy)= z^ys k /kl 
fc= 0 

The n columns of the Jacobian matrix of /, df/dx, are given by 

^ = ({/[(*, -■,{/[(*. Mil) (3-31) 

where we are again using the shorthand, for any vectors x, z 

(x,z) = (x,z, 0,...,0) 

and 6i is a unit column along coordinate X{. Thus n first-order calls of the vector-function algorithm 
produces its Jacobian matrix. 

Examples 

Scalar- and vector-form functions are useful as building blocks of higher level algorithms. 


Example 3.1. Spherical and Cylindrical Coordinates. Consider the back and forth transformation between 
Cartesian coordinates and cylindrical coordinates and their m time derivatives. From the standard relations among 
coordinates, 

/ r cosV>\ 


OCj* — 


r sin 


(3-32) 


\ * ) 
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and conversely. 


f r\ 

W 


( \R l +*?2 \ 

arctan(x r 2 , x r \) 




Xr 3 


(3-33) 


y 


It follows that for cylindrical coordinates the dynamic forms transform according to the algorithm XY -<CY: 

fSF m (r)* cos SF m (rl>)\ 

SF m (r) ★ sin SF m (ip) 

SF m (z) 


VF m (x r ) = 


\ 


(3-34) 


/ 


and conversely XY >- CY : 


( SF m (r)\ 
SF m (xP) 

V SF m (z) ) 


( [I/F m (x r )f + VF m {x r )l\ 1 / 2 \ 


arctan [FF m (a: r ) 2 , V F m (x r ) i ] 

VF m (x r ) 3 


(3-35) 


The conversion to spherical coordinates is similarly established. From the standard relations between coordinates, 
where 7 , xp are the latitude and longitude angles, respectively, 


x r = 


( p cos 7 COS 
p cos 7 sin ip 
\ -p sin 7 / 


(3-36) 


and, conversely. 




( \x r \ > 

V' 

— 

arctan(x r 2, x r i) 

\ 7 / 


\arctan(— x r 3,cos^x r i + sin^ r 2 ) ) 


(3-37) 


It follows that for spherical coordinates the dynamic forms transform according to the algorithm XP -<SP: 

( SF m (p) * cos 5F m ( 7 ) * cos SF m (xP) \ 


VF m (x r ) = 


SF m (p) * cos SF m ( 7 ) * sin SF m (xp ) 
—SF m (p) * sin 5 'F m ( 7 ) 


(3-38) 
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and, conversely, XPySP: 


SF m (p) \ 


SF m {tp) 

= 

VSF W (7)J 

\ 


\SF m (x r )\ 

arctan[VF m (x r ) 2 , VF m (x r )i] 


(3-39) 


Thus we have the following bidirectional link between the Cartesian (XY), spherical (SP) and 
cylindrical (CY) coordinates. 

SP <=> XY <=> CY 
Of course, now it is possible to concatenate: the algorithm 

CYySP = XYySP(XY^CY) (3-40) 

transforms cylindrical to spherical coordinates and their m derivatives by calling first the algorithm 
XY -<CY and then algorithm XY y SP. 


Example 3.2. Air Velocity. Suppose that the cylindrical coordinates of the wind field depend only on altitude, 
that is, the xy-magnitude, the xy-direction, and the vertical components are given by 


fv w \ 

\™ w J 


= w(h) 


Then the Cartesian coordinates of the wind evolve according to 


VF m (w r ) = XY -< CY{w[SF m (/i)]} 


If an aircraft flies along a trajectory VF m+1 (x r ), then its altitude evolves according to (for aircraft, by conven- 
tion, the z -axis points down) 

SF m (h) = -VF m (xrh 


and its air velocity 

Its airspeed v a , glidepath angle 7 “ 
given by 


VF m «) = VF m (x r) - VF m (w r ) 

, and heading angle ip a , as well as their time derivatives up to order m, are 


( SF m (v a ) \ 


SF Tn (ip a ) 


XYySP[VF m (v £)] 


\SF m ( 7 a ) ) 


Next, we consider rotations. 
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4 ROTATIONAL FORMS 


Rotational forms are useful for describing the attitude of a rigid body, its angular velocity, and its 
derivatives. 

Let a and 6 be two arbitrary (right-handed orthonormal) axis systems. The *a coordinates of b 
will be denoted by the matrix C ba . Thus, for i = 1, 2, 3, 

3 

= 53 Cba (h (4-1) 

3 = 1 

Consequently, the coordinates of a vector x transform as follows: 

x b = C ba x a (4-2) 


Since C^i, j) = 6j • fij, C ba is a direction cosine matrix of 6 relative to *a. The j th column of 

Cba gives the 6-coordinates of ay, the row, the a-coordinates of 6j. A direction cosine matrix is 
orthogonal so that its inverse is given by its transpose 

C ba = C ba= ^b (4-3) 

Since a and 6 are both right handed, the determinant 

detCft a = +1 (4-4) 

So, C ba always has an eigenvalue of +1. The corresponding eigenvector is given by axis (Cba) (provided 
axis^fa) ^ 0) where, for any 3 x 3 matrix C, the axis function 


axis(C ) = — 


( c 23 - c 32 \ 
C 31 - c 13 
\ci 2 - C21 / 


It may be noted that for any three-dimensional x. 


(4-5) 


axjs[5(x)] = x (4-6) 

and for any 3x3 matrix C, 

S'[axis(C)j = 1/2 (C - C T ) (4-7) 

We shall denote the set of all rotations as well as all 3 x 3 direction cosine matrices by 50(3). 

Suppose that C ba is a function of time. Then, since = /, it follows that Cb a C ba + 

CbaCl= 0. Therefore, Cb a C ba is skew symmetric, so let 


Cba Cfoa ~ S((jJbab) 


(4-8) 
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Consequently, 


u bab = axis(Cb a Cb a ) 


(4-9) 


and 

c ba = SMCto ( 4 - 10 ) 

The 3x1 column matrix u>b a b gives the b coordinates of the angular velocity of b relative to a . The 
last subscript in u ba b denotes the coordinate system. Thus, u>b aa — C a b^bab gives the a coordinates of 
the same velocity. On the other hand, u a b a = —u>i )aa are the *a coordinates of the angular velocity of 

a relative to b . Also from equation (4-2) 

Xb = Cba x a "T CbaXa — Ri.Ubab^^'ba'^Ci 4” C b aXa Riubabi^b 4” C b aXa 

Thus, we obtain the Coriolis derivative, which computes the derivatives of the b -coordinates of x from 
the derivatives of the "a -coordinates of x and the angular velocity of b relative to a : 

Xb = Si^bab^^b 4 " ^-'baP^o. ( 4 - 11 ) 


The matrix form corresponding to the direction cosine matrix Cb a locating b relative to a will be 
defined by that is, 

MF Z=(c ba ,C ba ,...,C t £ ) ) (4-12) 

Often angular velocity and its time derivatives are of interest. Hence, we introduce another type of 
dynamic forms, namely, rotational forms. The rotational form of order m for attitude C ba will be 
defined as the direction cosine matrix and angular velocity together with its time derivatives: 


RFB = < 


Cba m = Q 

l(Cfc.,VF< m - 1 >(w w ,)) m> 0 


(4-13) 


where attitude Cb a € 50(3) and oJbab € R 3 gives the 6 -coordinates of the angular velocity of b 
relative to a . 


Algebra 

The identity is given by 

RF™ = (1, 0) (4-14) 

Consider the transformation MF -< RF constructing the matrix form that corresponds to a given 
rotational form. Since Cba = s { u bab) c ba^ MF m (C ba ) is given iteratively by 

cL‘ +1) = £ cf s(“£r’ , ) c £ ) (4 - |5) 

i=0 


Conversely, it follows from 

u ba b = axis(C ba C £) 
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that the function MF y RF making the rotational form corresponding to a matrix form is given 
iteratively by 

"L +1) = E Cfaxisicg-^ci?) (4-16) 

t=0 

Note that C ab is the transpose of C ba . Thus, we have the pair, 

' RFfc = MFyRF(MFfc) 

< (4-17) 

^ MFfZ = MF^RF(RF^) 

Multiplication of rotational forms may be defined by 

RF% * RF£ = MFyRF[MF -< RF(RFJg) * MF^RF(RFfc)] (4-18) 

and inversion by 

(RF%)~ 1 = RF% = MF y RF{[MF^RF(RFj£)] T } (4-19) 

where the transpose acts on a matrix form and so it is already defined. Consider a sequence of rotations. 
Suppose that there are three coordinate systems, *a, b , and *c. The rotation of *c relative to *a is given 
in terms of rotations c relative to b and b relative to a by the product of forms, 

RF™ = RF%*RF% (4-20) 

It may be noted that the rotation and matrix forms behave algebraically just as their generating direction 

cosine matrix. 

A generalization of the Coriolis derivative (4-11) may be obtained as follows. Let VF m (x r ) be the 
vector form describing the motion of vector x to order m with respect to r . It and its time derivatives 

up to order m with respect to another set of axes, b , which is rotating according to RFJ£, are given by 
VF m (x b ) = RFg * VF m (x r ) = MF^RF(RFg) * VF m (x r ) (4-21) 


Example 4.1. Angular Acceleration of Aircraft Stability Axes. As an application of the algorithms just 
developed, consider the following problem, which may occur in the synthesis of reference trajectories for an 
aircraft. Suppose that the velocity of the aircraft relative to the runway r (assumed to be inertially fixed) is 
given as a function of time, v r (t). Consider an axis system v aligned so that t7i is parallel to the velocity tJ and 
the total nongravitational force is in the v\-V 3 plane with negative projection on Find the direction cosine 
matrix C vr , locating v relative to r as well as the angular velocity u> vrv , and its two time derivatives and the 
force f v (in g’s) and its derivatives. We translate the problem in the usual way as follows ( g is the acceleration 
of gravity). 

Vl = v/\v\ 


f = g 1 v-r 3 
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V2 = V\X // |vi X f\ 


V 3 =V 1 X V2 



' vr 


V 


T 
2 r 


\vlJ 


Now, in order to get (uj vrv ,u; vrv ,u > vrv ), three time derivatives of these expressions must be taken. Without 
dynamic forms this task is difficult, but with dynamic forms it is easy. We simply replace any implicit or explicit 
variable of time in the above set of equations with the corresponding forms in r : 


VF 3 (v lr ) = VF 3 (v r )/\VF 3 (v r )\ 


VF 3 (f r ) = g~ l VF 3 (v r ) - VF 3 (6 3 ) 


VF 3 (v 2r ) = VF 3 (v lr ) x VF 3 (f r )/\VF 3 (v lr ) x VF 3 (f r )\ 


VF 3 (v 3r ) = VF 3 (v ir ) x VF 3 (v 2r) 


fVF 3 (v lr ) T \ 


mf 3 t = 


VF 3 (v 2r ) T 


\VF 3 (v 3r ) T J 


RF 3 r = RF ~< MF(MF 3 r ) 

The coding is done! The direction cosine matrix, the angular velocity and its two derivatives can now be read 
off RF 3 r (see appendix C). 


The example shows how dynamic forms may be used easily to obtain time derivatives to arbitrary 
order from complicated expressions involving functions of scalars, vectors, and matrices. The dynamic 
forms keep track of many variables and many detailed computations. In the next two sections, we 
consider Euler angle forms and Euler parameter forms and derive the appropriate transformations. 
Euler angles and parameters are discussed in reference 1. 
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Euler Angle Form 


In this section we develop the direct and inverse functions relating Euler angles and m of their 
time derivatives to the corresponding (direction cosine) matrix form. Suppose that is a function of 
time and that u bab = u, a constant unit vector. Then equation (4-10) becomes 


Cba = S(u)C ba 


(4-22) 


Consider this to be a differential equation with (7^(0) = I. It is linear with constant coefficients, so 
the solution at t — <f> is 

CbaW = e S( " W (4-23) 

But 

e S(u)4> _ j + + S 2 {u)<t> 2 /2 + . . . (4-24) 

which becomes, on repeated application of the identity S 3 (u) = —S(u), 


e S{u)4> _ / + s in0S'(u) + (1 - cos <p)S 2 (u) = cos <\>I + sin (f>S(u) + (1 — cos <t>)uu T 
Note that, since S(u)u = 0, 


= u 


that u is the Euler axis of C fa, that 


and that the trace 


axis(e^^ u ^) = sin0u 
tr(e s M*) = 1 + 2cos0 


(4-25) 
(4-26) 

(4-27) 
(4-28) 

According to the Euler theorem on rotations, any attitude Cfa. may be achieved (parameterized) by a 
single rotation from / about an axis u through an angle <j>. 

An elementary Euler rotation about axis 8 q for q— 1,2,3 through angle a is given by 

Eq{a) = = cos al + sin a S(S q ) + (1 — cosa)6 q 6j (4-29) 

or explicitly: 


(4-30) 



(' 

0 

° \ 

Ei(a) = e S ^ a = 

0 cos a 

sin a 


Vo - 

sin a 

cos a) 


/cos a 

0 

— sina\ 

E 2 {a ) = e s ^ a = 

0 

1 

0 


V sin a 

0 

cos a / 


(4-31) 
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cos a sin a O' 


E^(a) = e S ^ a = — sin a cos a 0 


0 1 , 


The nine possible columns are given by 


Eq(ot)8j — 


if j = q 


cos aSj + sin aS(b q )bj else 


The following index functions for the index set {1,2,3} will be used to simplify notation. 


For i ^ j 


ix{i,j) = i — (j mod 3) 


Hhj) = 


-1 if = 2 


(4-36) 


Then, for example, 

SiSJSj = h(i,j)6 HiJ) (4-37) 

The function v(i,j) returns the integer that is different from either i or j. The function h(i,j) returns 
the sign of the projection of S(S{)6j on It may be noted that 


h(j,i) = ~h(i,j) 


(4-38) 


Wth the aid of these index functions equation (4-33) may be changed to the following more convenient 
form 


r <5 j if 3 = q 

Eq(oc)6j = l 

{ cos aSj + h(q, j) sin ab^ q ^ else 
Therefore, if C = E q (a), the elements of C are given by 


(4-39) 


Cij = 


if 3 = q 


cos ab ij + h(q, j ) sin ab iv ^ q ^ else 


(4-40) 


where axes subscripts on the direction cosine matrix C have been dropped temporarily to simplify 
notation. It is easy to raise this algorithm to order m: 
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ALGORITHM: MF m {C ) = Fg[SF m (a)] 


do 1=1,3 
do j = 1,3 

MF m (C)ij = 


end do 


r SF m {6ij) 

k cos SF m {a)6 ij + h(q,j ) sin SF(a)6 iu{qj) 


if j = Q 


else 


end do 


Next consider sequences of elementary rotations. Let C be a composite of three elementary 
rotations, with q\ ^ <72 / <73, that is, an Euler sequence 

C = E qi (ai)E q2 (a2)Eq 3 (a^) (4-41) 

Note the usual reverse order. The first rotation in the sequence is £^3(03): it is about axis 93 through 
angle 0:3. If q\ = 93, we will refer to the Euler sequence as repeating; if q\ ^ 93, nonrepeating. The 
elements of C may be computed by the indicated matrix product in equation (4-41). See page 20 of 
reference 1 for an expanded view of equation (4-41) for the twelve possible sequences. 


Example 4.2. Normally in flight control the attitude of the aircraft body axes b relative to the runway axes r 
may be represented by the (1,2,3) Euler sequence; 

C br = E 1 (4>)E 2 (d)E 3 (xP) (4-42) 

Alternatively, aircraft attitude may be represented by a redundant sequence: 

Cbr = E 2 (a) Ej {/3)E 1 (a)E 1 {<t> v )E 2 (iv)E3{ipv) (4-43) 

The first two rotations define the wind axes w> in which w\ is along the relative-velocity vector 

v a = v — w, where v is the aircraft velocity, w is the wind velocity, and W 2 is horizontal, that is, i? 2T3 = 0. The 
roll through angle cj> v defines the aircraft stability axes v . The body axes are reached by an additional corrective 
roll a, sideslip /?, and angle of attack a . The aerodynamic forces and moments are typically given as functions 
of airspeed |u a |, a, f3, and other variables. 


Consider again equation (4-41). Since the algorithms for E q [SF m (a )] and ★ have already been 
constructed, the zero-order algorithm (4-41) may be raised to order m: 

MF m (C) = £„[SF m (a 1 )] * £, 2 [SF”*(a 2 )] .£, 3 (SF"> 3 )] (4-14) 

This algorithm will be denoted as 


where the Euler angle form 


MF m (C ) = MF -< AF[AF™(a)\ 


AF™(a) - 


VF m (a ) 
Q 


(4-45) 


(4-46) 


¥ 


m 
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where the three angles are represented by a column. 


/“ 1 \ 


a = 


a 2 


W 


(4-47) 


and where the rotation sequence is represented by a row, 

g = (91,92,93) 


(4-48) 


Conversely, consider the inverse problem of computing the Euler angles from the direction cosine 
matrix. Consider first the middle angle 02 - From the defining equation (4-41) it follows that 


£92 ( a 2) = Eqi(-<Xl)CE q3 (-a3) 

and, since 6 q \ and S q3 are eigenvectors of E qi and E q3 , respectively, that 

6q 1 E q2 (a 2 )6 q3 = Sj x [cos a 2 6 q3 + h(q 2 ,q 3 ) sin a 2 6^ q2m) ] = 6 qi C6 q3 
If 91 = 93, *en, since = 1 and S qx K(q 2 ,q 3 ) = °* let 

x 2 = Cqiqx 


2/2 = \ A “ x 2 


If 91 # 93, then, since 6^6 q3 = 0 and $ qi S v ( q2m ) = 1. let 

2/2 = ^(92, 93)^193 
x 2 = y{l-\ 


v 2 

2/2 


Then, in either case, 

The range of a 2 is given by 


R 2 = 


a 2 = arctan (7/2,22) 


{0 < a 2 < 7r } , 


if 9i = 93 


{-n/2 <a 2 < tt/ 2}, if q\ ^ 93 
The solution corresponding to the negative square root is 

-a 2 if 9l = 93 
-q 2 + 7 r if 91 7^ 93 


„ * 

a 2 = 


(4-49) 


(4-50) 


(4-51) 


(4-52) 


(4-53) 


(4-54) 


(4-55) 
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The range ol\ is given by 


m = < 


r {-7r<a^<0}, if 91 =93 

< 

< {tt/2 < c *2 < k ] U {— 7T < 0-2 < — 7r/2}, if gi ^ 93 
Conditions on the boundary separating regions R 2 and iftjj are known in practice as gimbal lock 


(4-56) 


Next consider angle a\. From the defining equation (4-41) it follows that 

E q2 (a 2 )E q3 (a 3 ) = Eq x {a\)C 

so that 

^92 E( i2 ( a 2 ) Eq^ (a,)6 Q3 = 0 = 8q 2 Eq x (ai)C6g 3 


Equation (4-39) may be used to simplify the premultiplier of the last term: 

0 = 6g 2 Ej x (ai)C6 q3 = cosait/i - sinaixi 


where 


and 


X 1 92 )Ci/(<ji ,<72)93 

, 2/1 = a (?) ^9293 


f 1 


<7(9) = 


[ ^^(^<72)^73 


if 9l = 93 
if 91 + 93 


It may be noted that, for 91 7^ 93, 

<7(9) = h(qi, 93) = —h(q 2 , 93) = -h{qi,q 2 ) 


(4-57) 

(4-58) 

(4-59) 

(4-60) 

(4-61) 

(4-62) 


If |xi | + |t/i | = 0, then we have a singular case (gimbal lock) where the middle rotation makes the outer 
rotations equivalent: only the sum Q\ + 03 can be computed from C. As noted previously, gimbal 
lock occurs at ct 2 = ±n/2 for a nonrepeating sequence, and at a 2 — 0, it for a repeating sequence. 
Away from gimbal lock, there are two solutions depending on the choice of sign for components in 
equation (4-59). One solution is given by 


a\ = arctan(yi,xi) (4-63) 

and the other solution is given by 

ot\ = arctan(— yi, — xi) = ai + tt (4-64) 

The first choice corresponds to ot 2 , whose range is R 2 . Thus, if q\ ^ 93 and a = 0 so that C = I, then 

' xi = h 2 (q\,q 2 )C q3q3 = 1 
, y\ = ^(q)c q2q3 = 0 




M 
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so ax = 0 and not 7r. On the other hand, if q\ = <73 and a = (0, 7 t/2,0) so that C = E q2 (n/2), then 
(using eq. (4-62)) 

' xi = -h(qi,q 2 )6„( qi q2) E q2 (n/2)6 q3 = h(qi,q 2 ) 2 = 1 

i 

< y\ = ^92,93 = o 

Therefore, in either case we get back a\ and not a\, which corresponds to ajjj. Similarly, for 03: 


Let 


rp 

8 qi CE qz (-OLz)8 q2 = 0 = 2/3 cos 03 - x 3 sin a 3 


x 3 = -<r(QMQ2,q3)C qiMq2q3) 
{ 3/3 = °{<l)C q 192 


Then 

The other choice of sign leads to 


03 = arctan(r/ 3 ,x 3 ) 

<*3 = arctan(— 2/3, -x 3 ) = a 3 + n 


(4-65) 


(4-66) 

(4-67) 

(4-68) 


The companion to a 2 is a 3 , and the companion to is 03. Thus we have the following algorithm 
(ref. 27) for extracting Euler angles from direction cosine matrices: 


ALGORITHM: (C,q)y(dt) 
if qi = 93, then 

x 2 ~ ^9191 
y 2 = y/i -x| 
else 

V2 = 93) ^*9193 

x 2 = \l^-v 2 

end if 

XI = -<r(q)h(qi,q 2 )C „ {qi , q2)q3 
V\ = 0 ’( 9 )C , 9293 

x 3 = -<T{q)h{q2, 93 )^ 919 ^ 2 , 93 ) 
y3 = ^(9)<^9192 
ax = arctan(yi,xi) 

«2 = arctan(?/2,X2) 

e*3 = arctan(2/ 3 ,x 3 ) 


Example 4.3. Let q = (1,2, 3). (All twelve cases are given in table 2.1 in ref. 1.) Then 


C = Ei(ax)E 2 (a 2 )E^(a^) 


C2C3 C2S3 — «2 \ 

— C1S3 -I- siS 2 c 3 c i c 3 + s\s 2 s% s\c 2 
\ S1S3 + C1S2C3 — S1C3 + cis 2 s$ c\c 2 
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where we abbreviate cosai = c\, etc. Following the algorithm, h{ 2 , 3 ) = — 1 ; so 3/2 = — C13 = S2 and 

x 2 = yr^sf. 

Next, a(q) = + 1 , h(l, 2 ) - — 1 , and v{\, 2 ) — 3 ; so x\ — C33 = C1C2, and 3/1 = C23 = sic 2 . 

Finally, v( 2 , 3 ) — 1 ; so x 3 = — = C2C3 and 3/3 = C12 = C2S3. 

Tims 

ai = arctan(siC2,ciC2) 

l «2 = arctan(s2. Nl) 

. 03 = arctan(c 2 S 3 i c 2 C 3 ) 

If one begins with 0 so that 

( — 7r < 0i < n 
< — 7 r /2 < 02 < 7 r/ 2 
. — 7T < 03 <TT 

then the combined computation 

0 yCya 

will produce an a = 0 . For 0 2 in the other sector Rt,, the combined computation will produce an a = 0 *, so 
that the alternate solution a* = 0 . 


The general algorithm Cya may be raised to order m as follows: 

ALGORITHM: AF™{<x) = M F y AF q [M F m (C)] 

if q\ = 93, then 

SF m ( xj) = MF m (C) nqi 
SF™(y 2 ) = y'(l) - [SF™(i 2 )] 2 

s 1 s 0 

SF m (y 2 ) = h(g2,Qz)MF m (C)q l q 3 

SF™(x 2 ) = v'd) - [S'F>»(!/2 )] 2 

end if 

= -o-(g)%i,g 2 )MF™(C) v(9l ^ 

5F"»( yi ) = <T{q)MF m {C) q2qz 
SF m (x 3 ) = 

SF™(t/ 3 )=o^)MF m (C% ig2 
5F m (ai) = arctan[5F m (t/i), 5F m (xi)] 

SF m (ct 2 ) = arctan[5F m (3/2)1 SF m (x 2 )j 
SF" 1 ^) = arctan[5F m (t/ 3 ), SF m (x 3 )j 

We have just developed an algorithm for extracting Euler angles and their m time derivatives from 
the direction cosine matrix and its m time derivatives for any Euler sequence. In fact we have also 
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linked the Euler angle forms ( AF ) to the rotational forms (RF): 

AF q <*=> MF RF 

So, for example, the algorithm RF y AF q extracting the Euler angle form from a rotation form may 
be executed by first transforming the rotational form RFj£ into the matrix form MF m (C br )\ then the 
matrix form is transformed into the Euler angle form AF 9 m (%), namely, 

AF™^) = MFy AF q [MF -c RF(RF^)} (4-69) 

and, conversely, the algorithm for computing the rotational form from the Euler angle form is 

RF£ = MFy RF{MF -< AF^AF™^)]} (4-70) 

As a reminder, we note that angular velocity, angular acceleration, and higher time derivatives, which 
are of particular interest in practice, are contained in the rotational form 

*fS!-(Cir, VF™~'M) (4-71) 

so that, if we have computed RF 2J?, then we have also computed and m — 1 of its time derivatives. 


Example 4.4. Consider again example 4.1, where the aircraft trajectory represented by the vector form V F 4 (ty) 
was transformed into the rotational form RF$ r , which describes the rotation of the stability axes v relative to 
the runway axes *r. Suppose that, following convention, we represent the attitude of the stability axes relative 
to the runway as 

Cyr = F\ {4 > v)F 2 irtv)F‘i (^v) 


Then the roll, flight path, heading angles 


/ 4>v \ 


Otyr — 


and their time derivatives up to order 3 are given by 


lv 

\^v ) 


^(12,3) Kr) = MFyAF {w) [MF<RF{RFl r )\ 
Furthermore, if the body attitude Cf, r = C^ v C vr and, following convention, 


C bv = E 2 (a)E 3 (-f3)E 1 (p) 


then, for the already computed body rotation the evolution of the angle of attack, the (negative) sideslip 

angle, and the relative roll angle p 

f a \ 


a bv — 


-V 


\ p ) 
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and their time derivatives are given simply by 

^?2,3,l)Kv) = MFyAF {2 j A) [MF<RF(RFl*RF? v )) 
where all operations have been already defined. A numerical example is given in appendix C. 

Example 4.5. Suppose that a sequence q and Euler angles a br are given and that the Jacobian J q (a br ) relating 
angular velocity c u brb to the Euler angle rates a br is desired so that 

^ brb ~ Jqfabr^^brb ( 4 “ 72 ) 

Then the three columns of J q are given by three calls to RF -<AF q : 

Jq( a br) = ( {RF^AF q [(a br , 8 i)]} 4 , {RF -<AF q [(a br ,62)]}4> {RF -<AF q [(a br , 63)]}4 ) ( 4 - 73 ) 

where, as before, ( x , y) = (x, y, 0, . . . , 0). 

Example 4.6. Suppose a satellite is required to move so that the Euler angles for a sequence q are a given 
function of time, a br = f(t). In order to check whether the maneuver is executable, we wish to check the 
required angular acceleration and acceleration rate. Let 

VF 3 (a br ) = ,/(0,/ (3) (<)] 

The required angular acceleration and acceleration rate are given by the third and fourth items in 
RF br = (CbriUbrbiUbrb^brbi • • 0 = MF)~RF{MF -< AF q [AF% (a br )]} 

Note that the Jacobian is not needed for this computation. 

Example 4.7. Still another example is provided by patching two Euler angle parameterizations. Let one set of 
Euler angles a br be given in the nonrepeating sequence q = (1, 2, 3) and another set (3 br be given in the repeating 
sequence p = (3, 2,3). Both represent the same motion of the body RFj. [j?. The two coordinate systems (angles 
and m time derivatives) are related by the back-and-forth maps, 

AF™({3 br ) = MFyAF p {MF<AF q {AF™(a br )]} 

and, conversely, 

AF™(a br ) = MFyAF q {MF^AF p [AFp(P br )]} 

Such changes in coordinates are useful when the trajectory passes near gimbal lock in one set of coordinates. 


Euler Parameter Form 


Often in attitude control it is desirable to express attitude error in terms of Euler parameters, 

e = sin(0/2)u 


(4-74) 


[ r? = cos(0/2) 

where the unit column u is the Euler axis of rotation and cf) is the angle. The direction cosine matrix is 
given in terms of Euler parameters by 


C = I + 2 r]S{e) + 2 S 2 (e) 


(4-75) 
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and, conversely. 


r v = J[tr(C) + 1] 1 / 2 

< (4-76) 

k e = \axis{C)/r] 

The constraint is e 2 + rj 1 = 1. A singularity exists at 0 = 7r since there are two solutions ±u. The 
Euler parameter form is defined by the Euler parameters and their m time derivatives. 


ppm 


fVF m (e ) 
\SP m M 


(4-77) 


The relation between MFj£ and PFj£ is obtained simply by rewriting the defining equations in 
terms of dynamic forms. That is, the algorithm for 

MF$ = MF -< PF(PFfr) (4-78) 

is obtained simply by converting equation (4-75) to dynamic forms: 

MF m (C ) = MF m (I) + 2SF m {ri) * 5[KF m (e)] + 2S[VF m {t )] ★ S[^F m (e)] (4-79) 

and, conversely, the inverse algorithm 

PFfr = MFy PF(MFj%!:) (4-80) 


is the translation of equation (4-76): 

r SF m (ri) = tr[MF m {C )] + SF m ( l)} 1 / 2 

< (4-81) 

. VF m (e) = \axis[MF m {C)}/SF m (i}) 

The angle forms ( AF q ), rotational forms (RF), and parameter forms (FF) are now linked to the 
(direction cosine) matrix forms (MF) as follows: 

AF q & MF <=> RF 

t 

PF 

In other words, the algorithm for converting Euler angle form to Euler parameters form is given by 

PF{£ = MF >- PF{MF -< AF q [AF™ )] } (4-82) 

and similarly for any two representations of rotational motion shown in the diagram. 
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Attitude Servo 


Consider a simplified attitude control problem. Suppose that the object to be controlled (the plant) 
is spherically symmetric and that it is described by the following set of equations: 



( Cb r 



djdt 

^ brb 

— 

z + V 


\ Z J 


\ u / 


(4-83) 


where and are, respectively, attitude and angular velocity of the body b relative to a reference 
r . There is an intervening integral between the angular acceleration and the control u. In addition, 
there is a disturbance v. A block diagram of this plant is shown in figure 4.1. Note that the plant is 
of the pure feedback type. The objective is to devise a control scheme for tracking time-varying servo 
input Ccr(t)- The servo design shown in figure 2.8 may be specialized to attitude control as follows: 

Step 1. Command. First consider the command. Construct the corresponding third-order matrix 
form MF^. from the known function of time for the servo input Ccr(t ) for t > 0. For example, suppose 
that the command generator produces Euler angles and derivatives up to order three, AFg(acr). Then 

MFl = MF^AF q [AF^{aar)} 

Alternatively, M F^ r may have been computed from a translational trajectory as in example 1 where the 
dynamic form of the stability axes was computed. 


Step 2. Output. Next, from the plant state estimate (assumed to be available at all times t > 0), 


( ^br \ 


x = 


^brb 


V £ ) 


(4-84) 


and an estimate of the disturbance and its rate V F 1 (u) with c = z + £>, construct the plant output 
as follows: First compute an estimate of the plant rotation form. 


RF br = {Cbri 


& brbi ^ brb ) 


(4-85) 



v 

Figure 4. 1 . Plant block diagram. 


T 
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Then transform it into matrix form, 


M F l = (c*r, <4 C br )=MF*RF(RFl r ) (4-86) 


Step 3. Error. Now define the tracking error by the relative transformation, as in reference 28, 


MFl = MFl * MF? C 


(4-87) 


The goal is to drive the error close to ( /, 0, 0 ) and keep it there by manipulating the last element of 
( 3 ) 

RFfo , namely which is directly accessible for control. The control law for the error may be 
based on various parameterizations, such as Euler angles and Euler parameters. Following reference 29, 
we choose the Euler parameters (see eq. (4-81)), 

/e( 0), €<»), eW\ 

PFl=\ I = MF y PF(MFu) (4-88) 

W°>, W n<v) 

and view the vector part VF 2 (e) as the error state equation, where u e is written in place of e^ 3 ) to 
emphasize that e^ 3 ) is a control variable 



/ e (0) \ 


e (1) \ 

d/dt 

e(D 

— 

e< 2 > 




\ u e ) 


(4-89) 


The equation is a set of three strings each three integrators long (i.e., a Brunowsky form with Kronecker 
indices (3, 3, 3)), as shown in figure 4.2, where each e® is three dimensional. The control u € is as yet 
undefined. We design an asymptotically stable regulator law g in the usual way (additional dynamics 
such as integral feedback could be added to the control law if desired). 


U e 


S ( £ < ( V 1 >,e< 2 >) 


(4-90) 


Step 4. Plant Control. In order to implement this control law, u ( must be transformed back into 
the plant input u. Thanks to the regulator law, at this point we have raised the order m of VF m {e) 
from two to three: 

VF 3 (e) = [VF 2 {e),u e ] (4-91) 

and we impose the Euler parameter constraint, 

SF 3 (i]) = [SF 3 (1) - VF 3 (e) ■ VF 3 {e)} 1 / 2 (4-92) 



Figure 4.2. Brunowsky form for error dynamics. 
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thereby raising the order of PF\ ™ 

to three. Then we compute the new 



MFf c = MF -< PF(PFb c ) 

(4-93) 

- — - 77 ) 

and change coordinates to MF^ r : 

MFi = MFl*MFl 

(4-94) 

Finally, we compute the expanded 

rotation form. 



RFir = RF<MF{MF\ r ) 

(4-95) 

. — . 3 r. 

The last column of RFb r is 

Hence, the new value of the plant input is given by 



U = &brb ~ V 

(4-96) 


Thus, we have designed a control law for the attitude servo, which will track the variable input attitude 
G'bc(t) provided that the total error angle (j) < tt (eq. (4-74)) and that the control constraints are not 
violated. 


The block diagram of the resulting servo is shown in figure 4.3. 

The input may be given in terms of Euler angles with a fixed sequence, and the sequence may 
also be switched. In addition, Euler parameters and direction cosines may be used in any combination 
during operation. In all cases the input is transformed into the matrix form MF ^.. 



GENERATOR REGULATOR TRANSFORM PLANT 

Figure 4.3. Attitude servo block diagram. 
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If the moment-of-inertia matrix J b of the plant is not spherically symmetric, the gyroscopic terms 
may be viewed as a disturbance, 

^ J l) ^ (uJbrb) Jb^brb 

and, as before. 


U> 


(k+ 1) _ Ak) 

brb ~ z 


+ 


*(*> 


so VF^iu) can be easily computed by the techniques already described. 


The servo control algorithm is summarized in table 4-1, where n z is the number of integrators 
between u and which is in this case one. The algorithm is easily generalized to cases with n z > 1 
by making the order a variable, that is, by taking the highest order to be not 3 but 2 + n z . A numerical 
example of a scanning, spinning satellite is given in appendix D. 


Table 4-1. Servo control algorithm n z = 1 


Purpose 

Algorithm 

Given input at t = t n 

AFq (otcr) 

Change input to matrix form 

MF% = MF^AFq[AFq(ctcr)] 

Estimated plant state at t = t n 

(C br , Ufyrh ) 

Estimated disturbance at t = t n 

VF l (£>) = (£>,£) 

Estimated plant rotation form 

RF lyp = ( Cjryp , Wbrbl Z + V ) 

MfI = MF <RF{RF'\ r ) 

Change output to matrix form 

Change coordinates 

MF? c = M~Fl*MF? c 
PFl = MF^PF(MFl) 

Change to Euler parameters 

Error state 

(e, €, e) = VF 2 ( ( ) 

Raise order with regulator law 

C W=g(t, i, i) 

Error state and control 

VF 3 (e)=( e , i, e, e ( 3 >) 

Enforce constraint 

SF 3 (v) = [SF 3 { 1) - TF 3 (e) • VF 3 (e)] 1 / 2 

Change error to matrix form 

MFh — MF <PF(PF 3 c ) 

mPI = mf 3 c *mfI 

RF^ = MF y RF(MF br ) 

Change coordinates 

Change to rotation form 

Control at t — t n+ \ 

u = ( RFlh ~ i' 
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5 CONCLUSION 


The formalism of dynamic forms has been presented. The formalism may be used to translate 
effectively and systematically zero-order algorithms into m-order algorithms. In this way, for example, 
a large subroutine that computes the total force and moment acting on an aircraft for a given input, such 
as aircraft state, controls, and wind, may be routinely converted into the corresponding subroutine that 
computes the forces, moments, and time derivatives up to any order for given input and time derivatives 
to that order. Similarly, derivatives can easily be passed through elaborate coordinate changes. It was 
shown how this capability may be used to organize and simplify the design of control systems. Whereas 
many examples were provided to illustrate the application to automatic control, the emphasis was on 
the formalism of dynamic forms. Current research is concerned with specific aircraft applications and 
with the interpretation of control problems, in general, in terms of dynamic forms. 
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APPENDIX A 

NOE NUMERICAL EXAMPLE 


A specific numerical case of the nap-of-the-Earth (NOE) maneuver discussed in example 2.2 is 
given in this appendix. Let the height h and width a in feet of the Gaussian cap 

h = hmaxe-W^-W^ 2 


above a flat terrain be 


/ hmax N 


50 \ 


= 

50 

\ <Ty ) 


U00 / 


and let the horizontal path be as shown in figure A-l. The peak is at 50 ft; the level curves are 10 ft 
apart. The helicopter slides along the cap from the third quadrant, scoots around the peak, and leaves 



Figure A-l. Plan view of the flightpath and the Gaussian cap. 
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the area. The maneuver consists of four control points linked by polynomial segments each lasting 
4 sec. 

The scalar form of arc length SF^(s) is evolving according to figure A-2. The arc length s shown 
in the first panel grows from 0 to about 400 ft. The speed s shown as the solid line in the second panel 
is initially 50 ft/sec. During the first 4 sec it is reduced to 20 ft/sec, where it is held during the turn. 
Finally the helicopter accelerates at a constant 10 ft/sec 2 . The acceleration $ is shown dotted in the 
second panel. The first and second derivatives of acceleration are shown in the third panel, and 
the fifth time derivative of the arc length is shown in the bottom panel of the figure. The higher 
derivatives may be of interest for the verification that neither the pitch control nor the pitch-control 
rate exceeds its limit. Thus, if the approximation is made that longitudinal acceleration is controlled by 
the pitch angle, say 9, namely s = ~g0, then 9 = —s^/g determines the pitch-control requirement, 
whereas 9 ^ = — s^'/g determines the pitch-control rate. 

The scalar form of the heading angle 5F 3 (V0 is shown in figure A-3. Initially the heading is 
-30 deg, and it is held constant for the first 4 sec. Then follow two successive 90-deg turns to the 
right. Finally, the helicopter exits at 150 deg. The angular rate is bounded by 1 rad/sec, and angular 
acceleration, by 1 rad/sec 2 . The time rate of angular acceleration is shown in the last panel. This 
derivative is useful for the verification that the yaw-control rate is not excessive during the maneuver. 

The shape of the Gaussian cap and the evolution of the arc length and the heading angle are 
considered in this example to be independent inputs. The helicopter is constrained to stay on the 
Gaussian cap while executing the maneuver. Consequently the altitude becomes a dependent variable, 
whose evolution SF 5 (h) is shown in figure A-4. The vertical acceleration, which is bounded by 
16 ft/sec 2 , has a major influence on the required thrust. The higher derivatives indicate the dynamic 
requirement from the thrust generator. 

Figure A-5 shows the (scaled) energy flow 5F 2 (e). The total kinetic and potential energy e, the 
power e, and the time rate of power e are all shown. 

This example illustrates the ease with which a variety of time derivatives of transformed functions 
of time can be computed by means of dynamic forms. It should be noted that, since we are discussing 
command-trajectory generation, the differentiation of noisy data is not an issue. 
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Figure A-4. Evolution of the resulting altitude scalar form SF 5 (h). 
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Figure A-5. Evolution of the resulting energy scalar form 5F 2 (e). 
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APPENDIX B 


AUTOMATIC CONTROL EXAMPLE 


This appendix contains an application of the control structure shown in figure 2.8 and the corre- 
sponding control algorithm Cj, which is repeated below for convenience. 


ALGORITHM: U p = C f[x p , SF n (x i c )] 

SF n -\x lp ) = T lf (x p ,0) 

SF n ~ l (e) = SF n ~ 1 (xi p ) - SF n ~ l (x lc ) 
e( n ) = k[SF n ~ 1 (e)J 
SF n (xi p ) = SF n (xi c ) + SF n (e) 
[SF°(u p ),x p ] =T- f l [SF”(x lp )} 
u p = [SF°(u p )) 0 


Suppose that the 
follows: 


plant is described by a four-dimensional state and that the state equation is as 

( x 2p \ 



sin(0.2xi p )x2p + (2 + cosi)x3 p 

X4 p 


\ U P 


) 


This system is the same one used in example 2.3, where the algorithm was constructed, and in 
example 2.4, where the algorithm T^} was constructed. The only item remaining to be specified is the 
regulator. We choose a simple linear, constant-gain regulator: 


e = — k\e — 


and place one critically damped pole pair at 0.5 rad/sec and the other at 1 rad/sec. Of course limiters 
and additional dynamics such as integrators could have been used. Now Cf is completely specified. 


Next consider the inputs. Suppose that the input to be tracked consists of four segments, each 
lasting 4 sec and taking the system from x\ c = 0 at the beginning of the maneuver to 


/64\ 


xi c (16) = 


0 

0 


v 0 y 


at the end. The time history of the input SF^{x\ c ) is shown as solid lines in figure B-l. The system is 
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time (sec) 

Figure B-l. Evolution of input SF^(x\ c ) and response SF^(x \ p ). 
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. (2) (2) 
commanded to reach steady acceleration x\ c = 2 at the first waypoint at t = 4, then to reach x lc = 0 

at the second waypoint at t = 8 , where the speed x^ = 8 ; the boundary condition at the third waypoint 
( 2 ) 

is a steady x\ c = — 2 ; finally, the boundary condition at the fourth waypoint at t = 16 is the desired 
steady position x\ c — 64. The response of the plant SF‘ i (xi p ) with initial condition x p = 0 is shown 
in the figure by dotted lines. The higher plant state coordinates (x^ p , x^ p ) and the control u p are shown 
in figure B-2. Also shown in the figure are the ’’damping” term q\ = asin( 6 xi p )x 2 p as a solid line and 
the “effectiveness” of x^ p , namely q<i = [2 + cos(cf)], as a dotted line. The regulator error SF 4 (e) is 
shown in figure B-3. It may be noted that tracking is good despite considerable activity in q\ and q 2 - 
The region near t = 9 is especially disruptive because <71 is maximum while the effectiveness of x^ p , 
namely q 2 , is near minimum. 


The response to initial offset. 


/ 10 \ 


x p(°) = 


0 

0 


\ 0 / 


is shown next for the same input. Figure B-4, as figure B-l, shows both the reference and the response. 
The plant tracks the input asymptotically. In addition, as shown in figure B-5, the error behaves as an 
autonomous response of constant linear dynamics with the assigned poles. The sampling rate was 100, 
and the integration step size, 0.01. 


This simple example illustrates the application of dynamic forms to control. It should be clear that 
the same control algorithm Cf can be easily extended to practical cases in which, as noted previously, 
the implementation of the zero-order function / requires more than 4,000 lines of FORTRAN. An 
application of Cf to attitude control is given in appendix D. 
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Figure B-2. System response in terms of natural coordinates. 
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Figure B-3. Error response SF 4 (e) for x p (0) = 0. 
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Figure B-4. Evolution of input SF^(x i c ) and response SF 4 (x \ p ) for xi p (0) = 10. 
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APPENDIX C 


STABILITY AXES NUMERICAL EXAMPLE 


For a specific numerical case of example 4.1, we consider the Euler angles of the stability-axes 
system. Let 

/ cos 7 cos ip \ 


V r = VU r — V 


cos 7 sin ip 


\ -sin 7 / 

and let the speed v, the heading angle ip, and the glidepath angle 7 evolve according to 


/SF 4 (v)\ 


f V,v,.. 

. , ^ 


/100 

10 

-15 

20 

5 ^ 

SF 4 (iP) 

= 



= 

1.0 

0.3 

0.4 

-0.5 

-0.6 

SF 4 ( 7) J 


7.7. 

•, 7 ^/ 


V0.3 

0.1 

-0.2 

-0.1 

0 / 


These conditions may be at time t along the reference trajectory, which meets some given boundary 
conditions at the next waypoint. Then the evolution of the corresponding sines and cosines is given by 


/cos SF 4 (ip) \ 


COS Ip, (cos 'lp)( l \ . 

. . , (cOS^)( 4 ) \ 


/0.54 

-0.25 

-0.9 

0.25 

0.76 \ 

sin 5 F 4 (- 0 ) 


sinxp, (sin^)^ 1 ), . . 

, . , (sin^)( 4 ) 


0.84 

0.16 

0.14 

-0.59 

-0.33 

cos SF 4 ( 7 ) 


cos 7, (cos 7 )^, . . 

, . , (cos 7 )( 4 ) 


0.96 

-0.03 

0.05 

0.09 

-0.08 

sin 5F 4 ( 7 ) ) 


\ sin 7 , (sin 7 )^ 1 \ . . 

.,(sin 7 )( 4 ^ / 


\0.30 

0.10 

-0.19 

-0.08 

- 0 . 01 / 


Consequently, the path tangent vector u r evolves according to 



cos 5F 4 (7) ★ cos SF 4 (tp) \ 


l 0.52 

-0.26 

-0.33 

0.28 

0.45 \ 

VF 4 (Ur) = 

cos SF 4 ( 7) * sin SF 4 (xp) 


0.80 

0.13 

0.17 

-0.45 

-0.22 


— sin 5 F 4 ( 7 ) j 


\ -0.30 

-0.10 

0.19 

0.08 

0.01 / 


the runway coordinates of the velocity and their four derivatives are 



( 51 

-20 

-45 

40 

67 \ 



VF 4 {v r ) = SF 4 (v ) ★ VF 4 (u r ) = 

80 

21 

7 

-32 

-41 


* 


\ —30 

-12 

22 

12 

-22/ 
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and the runway coordinates of the total generated force and its three derivatives are (in g’s) 



/ —0.64 

-1.41 

1.25 

2.09 \ 

VF 3 (fr) = g~ l VF 3 (v r ) - VF 3 (6 3 ) = 

0.65 

0.22 

-1.01 

-1.28 


\ — 1.39 

0.68 

0.38 

-0.69/ 


For the unit vector along v r 



0.52 

-0.26 

-0.33 

0.28 \ 

^ 3 (t»lr) = VF 3 (v r )/\VF 3 (v r )\ = 

0.80 

0.13 

0.17 

-0.48 


—0.30 

-0.10 

0.19 

0.08 / 


It may be noted that (happily), VF 3 (v i r ) = VF 3 (u r ), as expected. The normalized cross product. 



/— 0.60 0.40 

-0.17 

— 1.65\ 

VF 3 (v 2r ) = VF 3 (v lr ) x VF 3 (f r )/\VF 3 (v lr ) x VF 3 (f r )\ = 

0.59 

-0.23 

-0.30 

1.53 


0.55 

0.68 

-1.10 

0.66 / 


and, for third-axis vector, 



/ 0.61 

0.61 

-0.86 

0.64 \ 

^(%) = X VF 3 (v 2r ) = 

-0.11 

-0.27 

0.95 

0.04 


V 0.78 

-0.51 

-0.09 

2.31/ 


The matrix form, MF 3 r = ( C vr , C vr , C vr , C'ffl ), is just a rearrangement of the vector form VF 3 (v ir ): 




0.52 

0.80 

-0.30N 


/— 0.26 

0.13 

-0.10\ 


MF 3 r = 


-0.60 

0.59 

0.55 

1 

0.40 

-0.23 

0.68 

5 ■ • • 



0.61 

-0.11 

0.78 ) 


0.61 

-0.27 

-0.51/ 



The ^-coordinates of velocity v and its derivatives are given by the Coriolis product, 

/100 10 -15 20 \ 


VF 3 (v v ) = MFl * VF 3 (v r ) = 


0 0 
\ 0 0 


0 

0 


0 

o / 
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Thus, as expected, the activity is only in the top row, which is equal to SF 3 (v). The ^-coordinates of 
the force / and its derivatives are given by 


VF 3 (f v ) = MFl*VF 3 (f r ) = 


( 0.61 -0.37 0.42 0.08 \ 

0 0 0 0 
V —1.55 -0.21 -0.06 3.90/ 


It may be noted that there is no lateral component: f v moves only in the v\ - plane, as required. 

Next, consider the angular velocity of *v relative to *r and its derivatives. We transform the matrix 
form M Fy r into the rotational form RF^ r by means of the function MFyRF to obtain: 




/0.80 

-0.98 

-0.20 \ 


RFy r — [ C vri V F^(uj vrv ) ] — 

Cvr i 

0.25 

0.21 

0.04 




\0.18 

0.20 

-0.37/ 



The direction cosine matrix C vr is the first matrix in MFy r , above, and 

2 

V F {uiyrv) = ( Wyrvi ^vrv > &vri j ) 


contains the angular velocity and its two time derivatives. The two derivatives of angular velocity are 
used to determine whether moment controls and control rates are within their limits. 


Finally, we obtain the Euler angles and their time derivatives for the sequence q — (1, 2, 3). 


AF% = RFy AF q (RFy r ) = 


/0.61 

0.89 

-0.83 

-0.33 \ 

0.30 

0.10 

-0.20 

-0.10 

1.00 

0.30 

0.40 

-0.50 

0 

1 

2 

3 / 


The first row corresponds to the evolution of the roll angle 5F 3 (0). The next two rows give SF 3 ( 7) 
and 5F 3 (^), respectively, which are the same (as of course they should be) as the spherical coordinates 
given at the beginning of this numerical example. 
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APPENDIX D 


ATTITUDE SERVO NUMERICAL EXAMPLE 


A specific numerical case of the attitude servo system developed in section 4 is presented in 
this appendix. The control algorithm for one intervening integrator is summarized in table 4-1. The 
regulator-control law chosen for this example is a simple, spherically symmetric (scalar gains) linear 
law with limited position feedback, namely 

u t = -sai(20.0e (0) , 1) - 18.14eW - 7.83e (2) (D-l) 

The feedback gains place the three poles at (s = —5, C = 0.707, u n — 2). The position feedback 
saturates at 1/sec 3 for each axis. Thus, saturation occurs for |ej| = sin(0/2)|uj| = 0.05, which, for 
ui = 1, is approximately 6 deg. The overall block diagram of the servo system is shown in figure 4.3. 

Now consider the input, given by the command generator as an Euler angle form AFq(acr). Let 
the Euler sequence q = (1, 2, 3), so that 

Car = Ei{a cr i)E2{a cr 2)E^{a cr z) 

The maneuver consists of four segments, each lasting 4 sec, as shown in figures D-l and D-2. 
Figure D-l shows the commanded motion in the yaw-pitch plane. Initially all three commanded angles 
are zero. The initial segment takes pitch a ^2 and yaw to 50 and -90 deg, respectively, and, as 
shown in figure D-2, brings the roll rate a^i to 1 rad/sec. In the second segment the pitch and roll 
rates are held constant while the yaw scans from -90 to 90 deg. During segments 3 and 4 the command 
remains at \ — 1 rad/sec, cst-crl = 50 deg, and 3 = 90 deg, while the second component of the 
disturbance v is raised from 0 to 0.1 rad/sec 2 , as shown in the bottom panel in figure D-3. The other 
components v\ = 1/3 = 0. This figure also shows the command in terms of the angular velocity and 
two of its derivatives. 

The state estimate is assumed to be exact; the disturbance is estimated erroneously by u = 0. The 
plant, which is initially offset from the origin by a^ r = (15, 15, —15) deg, responds as shown in 

figure D-4 in terms of Euler angles and in figure D-5 in terms of the angular velocity and its derivatives. 

The regulator error and Control SF^(t) are shown in figure D-6. It may be noted that the position 
feedback saturates, thereby reducing the error at constant rate, as shown in the second panel in the 
figure. The initial-condition transient is essentially over during the first segment. The final steady 
offset in €2 counteracts the disturbance in the usual way. Figure D-7 shows the error in a phase plane, 
where constant-speed sections are clearly visible. Note also that the changes in the input AFq(a cr ) are 
practically decoupled from the error. Finally, figure D-8 shows the tracking in the yaw-pitch plane with 
initial offset and disturbance. The hangoff in terms of Euler angle error is not steady (see fig. D-4). 
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Figure D-5. Evolution of estimated body rotation form RF^.. 
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Figure D-6. Evolution of error Euler parameter form PF^. 
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